# Bibliographic Information & Notes ############################################
# Title: Taking the Cloth -- Social norms and elite cues increase support for... 
# ...masks among white Evangelical Americans
# By: Claire Adida, Christina Cottiero, Leonardo Falabella, Isabel Gotti, ...
# ...ShahBano Ijaz, Gregoire Phillips, and Michael Seese
# Published In: Journal of Experimental Political Science
# Publication Date: 
# DOI: 
# Preprint: 10.31235/osf.io/yt3e7
# Note: This code replicates findings in the Supplementary Information Doc...
# ...for manuscript replication, please see TtC_Manuscript_Replication.R...
# ...code for SI-A.3, SI-A.4, and SI-B available in TtC_Manuscript_Replication.R


# Packages and Data ############################################################
## Working Directory ###########################################################
rm(list = ls())
setwd("/Users/mikeseese/Documents/GitHub/Taking_the_Cloth_Replication") # Change WD


## Packages ####################################################################
library("tidyverse")
library("modelsummary")
library("cobalt")
library("AER")
library("mediation")
library("MASS")
library("ggeffects")
library("emmeans")


## Data ########################################################################
working <- readRDS(file = "working.rds")


# Preliminary Data Cleaning ####################################################
## Recode Religiosity Index ####################################################
## Corrects a data transfer error from Lucid
working$rel1x <- dplyr::recode(working$rel1,
                               `1` = 1,
                               `2` = 6,
                               `3` = 2,
                               `4` = 3,
                               `5` = 4,
                               `6` = 5)

working$rel2x <- dplyr::recode(working$rel2,
                               `1` = 1,
                               `2` = 6,
                               `3` = 2,
                               `4` = 3,
                               `5` = 4,
                               `6` = 5)

working$rel3x <- dplyr::recode(working$rel3,
                               `1` = 1,
                               `2` = 6,
                               `3` = 2,
                               `4` = 3,
                               `5` = 4,
                               `6` = 5)

working$rel.ix <- working$rel1x + working$rel2x + working$rel3x


# Balance ######################################################################
## Balance Tables ##############################################################
## Demographic Balance
dem <- working %>% dplyr::select(
  `Condition` = treat,
  `Age` = age,
  `Gender` = gender,
  `Education` = edu,
  `Party ID` = party.id,
  `Follows News` = news.follow,
  `News Source` = news.source,
  `Follows Covid` = covid.follow,
  `Covid News Source` = covid.source,
  `Covid Threat` = covid.threat,
  `Tested for Covid` = covid.tested,
  `Friends with Covid` = covid.friends)

datasummary_balance(~Condition,
                    data = dem,
                    fmt = '%.2f')

## Religious Mechanisms Balance
mech <- working %>% dplyr::select(
  `Condition` = treat,
  `Regularly Attends Church` = rel1x,
  `Spiritual Values` = rel2x,
  `Americans Should be More Religious` = rel3x,
  `Religiosity Index` = rel.ix,
  `Religiosity (Self-Reported)` = religiosity,
  `Religious Motivation - Insurance` = c1, 
  `Religious Motivation - Morality` = c2, 
  `Religious Motivation - Social Capital` = c3, 
  `Religious Motivation - Heuristics` = c4,
  `Religious Motivation - Tradition` = c5, 
  `Religious Motivation - Social Identity` = c6, 
  `Religion Offers Comfort` = ei1,
  `Takes Religious Approach to Life` = ei2,
  `Enjoys Social Aspects of Church` = ei3,
  `Extrinsic / Intrinsic Motivation Index` = ei.index)

datasummary_balance(~Condition,
                    data = mech,
                    fmt = '%.2f')


## Love Plot ###################################################################
## Data Cleaning
dem <- dem[complete.cases(dem), ]
dem2 <- subset(dem, select = -c(Condition))

bt <- bal.tab(dem2, treat = dem$Condition, m.threshold = .1, 
              multi.summary = TRUE)

## Clean up variable names in Excel and save as .csv
var.names(bt, 
          type = "df", 
          minimal = FALSE)

nn <- read.csv("loveplot_newnames.csv")

## Love Plot
love.plot(dem2, treat = dem$Condition, 
                thresholds = c(m = .1), 
                binary = "std",
                which.treat = .all, 
                var.names = nn,
                abs = FALSE) + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  labs(subtitle = "Standardized Mean Differences Across Comparison Groups") +
  theme(legend.position = "none")

## Cleaning Up
rm(list=setdiff(ls(), "working"))


# Average Treatment Effects: OLS Tables ########################################
## Data Cleaning ###############################################################
w2 <- working %>% 
  filter(behave > 4)

iv1 <- " ~ treat + gender + age + education + party.f + covid.threat + covid.tested + covid.friends"
iv2 <- " ~ treat + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends"


## OLS Models ##################################################################
## Baseline Models
a1 <- lm(therm ~ treat, data = working)
a2 <- lm(p1.wash ~ treat, data = working)
a3 <- lm(p2.dist ~ treat, data = working)
a4 <- lm(p3.mask ~ treat, data = working)
a5 <- lm(p4.pray ~ treat, data = working)
a6 <- lm(behave ~ treat, data = working)
a7 <- lm(clicked ~ treat, data = w2)

## Demographic Controls
b1 <- lm(as.formula(paste("therm", iv1)), data = working)
b2 <- lm(as.formula(paste("p1.wash", iv1)), data = working)
b3 <- lm(as.formula(paste("p2.dist", iv1)), data = working)
b4 <- lm(as.formula(paste("p3.mask", iv1)), data = working)
b5 <- lm(as.formula(paste("p4.pray", iv1)), data = working)
b6 <- lm(as.formula(paste("behave", iv1)), data = working)
b7 <- lm(as.formula(paste("clicked", iv1)), data = w2)

## All Controls
c1 <- lm(as.formula(paste("therm", iv2)), data = working)
c2 <- lm(as.formula(paste("p1.wash", iv2)), data = working)
c3 <- lm(as.formula(paste("p2.dist", iv2)), data = working)
c4 <- lm(as.formula(paste("p3.mask", iv2)), data = working)
c5 <- lm(as.formula(paste("p4.pray", iv2)), data = working)
c6 <- lm(as.formula(paste("behave", iv2)), data = working)
c7 <- lm(as.formula(paste("clicked", iv2)), data = w2)


## OLS Tables ##################################################################
## Group models by outcome variable
therm <- list(a1, b1, c1)
wash <- list(a2, b2, c2)
dist <- list(a3, b3, c3)
mask <- list(a4, b4, c4)
pray <- list(a5, b5, c5)
beh <- list(a6, b6, c6)
click <- list(a7, b7, c7)

## Coefficients
cm = c("(Intercept)" = "Intercept",
       "treatEndorsement" = "Endorsement",
       "treatNorms" = "Norms")

## Adding Rows to the Tables
row <- tribble(~x0, ~x1, ~x2, ~x3,
               "Controls", 
               "No", "Demographic", "All" )
attr(row, 'position') <- c(7)


## Tables by Outcome Variable
modelsummary(therm,
             title = "Thermometer",
             coef_map = cm,
             add_rows = row,
             stars = TRUE,
             gof_omit = 'IC|Log|Adj')

modelsummary(wash,
             title = "Hand Washing",
             coef_map = cm,
             add_rows = row,
             stars = TRUE,
             gof_omit = 'IC|Log|Adj')

modelsummary(dist,
             title = "Distancing",
             coef_map = cm,
             add_rows = row,
             stars = TRUE,
             gof_omit = 'IC|Log|Adj')

modelsummary(mask,
             title = "Mask Wearing",
             coef_map = cm,
             add_rows = row,
             stars = TRUE,
             gof_omit = 'IC|Log|Adj')

modelsummary(pray,
             title = "Prayer",
             coef_map = cm,
             add_rows = row,
             stars = TRUE,
             gof_omit = 'IC|Log|Adj')

modelsummary(beh,
             title = "Would Click",
             coef_map = cm,
             add_rows = row,
             stars = TRUE,
             gof_omit = 'IC|Log|Adj')

modelsummary(click,
             title = "Clicked",
             coef_map = cm,
             add_rows = row,
             stars = TRUE,
             gof_omit = 'IC|Log|Adj')

## Cleaning Up
rm(list=setdiff(ls(), "working"))


# Compliance Average Causal Effects ############################################
## Generate Compliance Variables & Subset Data #################################
## Compliance Variables
as.data.frame(colnames(working))

working2 <- working %>% 
  unite("knew", 13:15, na.rm = TRUE, remove = FALSE)

working2$knew <- as.numeric(working2$knew)

working2$knewit <- dplyr::recode(working2$knew,
                                 `3` = "Knew It",
                                 `2` = "Didn't Know It",
                                 `1` = "Don't Believe It",
                                 .default = NA_character_)

working2$knewit <- as.factor(working2$knewit)
table(working2$knewit)

working2$comply <- dplyr::recode(working2$knewit,
                                 "Don't Believe It" = "Non-Compliant",
                                 "Didn't Know It" = "Compliant",
                                 "Knew It" = "Compliant")

table(working2$treat)
table(working2$treat, working2$comply)

## Endorsement Dataset
e <- subset(working2, treat != "Norms")
e$treat.e <- factor(e$treat, 
                    levels = c("Control", "Endorsement"),
                    ordered = TRUE)

e$t.e[e$treat.e == "Control"] <- 0
e$t.e[e$treat.e == "Endorsement"] <- 1
e$tu.e[e$comply == "Non-Compliant"] <- 0
e$tu.e[e$comply == "Compliant"] <- 1
e$tu2.e <- e$tu.e
e$tu2.e[e$t.e == 0] <- 0
# Coding Control as all non-compliant
# Because you can't comply with a treatment you didn't get
# tu2.e is now the compliance variabe: = 1 if compliant | endorsement condition

e2 <- e %>% 
  filter(behave > 4)

## Norms Dataset
n <- subset(working2, treat != "Endorsement")
n$treat.n <- factor(n$treat, 
                    levels = c("Control", "Norms"),
                    ordered = TRUE)

n$t.n[n$treat.n == "Control"] <- 0
n$t.n[n$treat.n == "Norms"] <- 1
n$tu.n[n$comply == "Non-Compliant"] <- 0
n$tu.n[n$comply == "Compliant"] <- 1
n$tu2.n <- n$tu.n
n$tu2.n[n$t.n == 0] <- 0
# tu2.n is now the compliance variabe: = 1 if compliant | norms condition

n2 <- n %>% 
  filter(behave > 4)


## Endorsement IV Models #######################################################
## Baseline
e1.1 <- ivreg(therm ~ tu2.e, ~ t.e, data = e)
e1.2 <- ivreg(p1.wash ~ tu2.e, ~ t.e, data = e)
e1.3 <- ivreg(p2.dist ~ tu2.e, ~ t.e, data = e)
e1.4 <- ivreg(p3.mask ~ tu2.e, ~ t.e, data = e)
e1.5 <- ivreg(p4.pray ~ tu2.e, ~ t.e, data = e)
e1.6 <- ivreg(behave ~ tu2.e, ~ t.e, data = e)
e1.7 <- ivreg(clicked ~ tu2.e, ~ t.e, data = e2)


## Demographic Controls
e2.1 <- ivreg(therm ~ tu2.e, ~ t.e + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = e)
e2.2 <- ivreg(p1.wash ~ tu2.e, ~ t.e + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = e)
e2.3 <- ivreg(p2.dist ~ tu2.e, ~ t.e + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = e)
e2.4 <- ivreg(p3.mask ~ tu2.e, ~ t.e + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = e)
e2.5 <- ivreg(p4.pray ~ tu2.e, ~ t.e + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = e)
e2.6 <- ivreg(behave ~ tu2.e, ~ t.e + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = e)
e2.7 <- ivreg(clicked ~ tu2.e, ~ t.e + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = e2)

## All Controls
e3.1 <- ivreg(therm ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e)
e3.1 <- ivreg(therm ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e)
e3.2 <- ivreg(p1.wash ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e)
e3.3 <- ivreg(p2.dist ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e)
e3.4 <- ivreg(p3.mask ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e)
e3.5 <- ivreg(p4.pray ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e)
e3.6 <- ivreg(behave ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e)
e3.7 <- ivreg(clicked ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e2)


## Norms IV Models #############################################################
## Baseline
n1.1 <- ivreg(therm ~ tu2.n, ~ t.n, data = n)
n1.2 <- ivreg(p1.wash ~ tu2.n, ~ t.n, data = n)
n1.3 <- ivreg(p2.dist ~ tu2.n, ~ t.n, data = n)
n1.4 <- ivreg(p3.mask ~ tu2.n, ~ t.n, data = n)
n1.5 <- ivreg(p4.pray ~ tu2.n, ~ t.n, data = n)
n1.6 <- ivreg(behave ~ tu2.n, ~ t.n, data = n)
n1.7 <- ivreg(clicked ~ tu2.n, ~ t.n, data = n2)

## Demographic Controls
n2.1 <- ivreg(therm ~ tu2.n, ~ t.n + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = n)
n2.2 <- ivreg(p1.wash ~ tu2.n, ~ t.n + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = n)
n2.3 <- ivreg(p2.dist ~ tu2.n, ~ t.n + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = n)
n2.4 <- ivreg(p3.mask ~ tu2.n, ~ t.n + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = n)
n2.5 <- ivreg(p4.pray ~ tu2.n, ~ t.n + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = n)
n2.6 <- ivreg(behave ~ tu2.n, ~ t.n + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = n)
n2.7 <- ivreg(clicked ~ tu2.n, ~ t.n + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = n2)

## All Controls
n3.1 <- ivreg(therm ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n)
n3.1 <- ivreg(therm ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n)
n3.2 <- ivreg(p1.wash ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n)
n3.3 <- ivreg(p2.dist ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n)
n3.4 <- ivreg(p3.mask ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n)
n3.5 <- ivreg(p4.pray ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n)
n3.6 <- ivreg(behave ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n)
n3.7 <- ivreg(clicked ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n2)


## CACE Tables #################################################################
therm <- list(e1.1, e2.1, e3.1, n1.1, n2.1, n3.1)
wash <- list(e1.2, e2.2, e3.2, n1.2, n2.2, n3.2)
dist <- list(e1.3, e2.3, e3.3, n1.3, n2.3, n3.3)
mask <- list(e1.4, e2.4, e3.4, n1.4, n2.4, n3.4)
pray <- list(e1.5, e2.5, e3.5, n1.5, n2.5, n3.5)
beh <- list(e1.6, e2.6, e3.6, n1.6, n2.6, n3.6)
click <- list(e1.7, e2.7, e3.7, n1.7, n2.7, n3.7)

row <- tribble(~x0, ~x1, ~x2, ~x3, ~x4, ~x5, ~x6,
               "Controls", 
               "No", "Demographic", "All", 
               "No", "Demographic", "All" )
attr(row, 'position') <- c(7)

cm = c("(Intercept)" = "Intercept",
       "tu2.e" = "Endorsement (Compliant)",
       "tu2.n" = "Norms (Compliant)")

modelsummary(therm,
             title = "Thermometer",
             coef_map = cm, 
             add_rows = row,
             stars = TRUE)

modelsummary(wash,
             title = "Hand Washing",
             coef_map = cm, 
             add_rows = row,
             stars = TRUE)

modelsummary(dist,
             title = "Distancing",
             coef_map = cm, 
             add_rows = row,
             stars = TRUE)

modelsummary(mask,
             title = "Mask Wearing",
             coef_map = cm, 
             add_rows = row,
             stars = TRUE)

modelsummary(pray,
             title = "Prayer",
             coef_map = cm, 
             add_rows = row,
             stars = TRUE)

modelsummary(beh,
             title = "Would Click",
             coef_map = cm, 
             add_rows = row,
             stars = TRUE)

modelsummary(click,
             title = "Clicked",
             coef_map = cm, 
             add_rows = row,
             stars = TRUE)

## Cleaning Up
rm(list=setdiff(ls(), "working"))


# Mediation Analyses ###########################################################
x <- " ~ treat + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"

## Religiosity Mechanism #######################################################
working <- rename(working, "religiosity.raw" = "religiosity")

working$religiosity <- factor(working$religiosity.raw)
working$religiosity <- dplyr::recode(working$religiosity,
                                     "1" = "Anti-Religious",
                                     "2" = "Not at all Religious",
                                     "3" = "Slightly Religious",
                                     "4" = "Moderately Religious",
                                     "5" = "Very Religious")

working2 <- working %>% 
  filter(behave > 4)

## Mediator
h4.med <- polr(as.formula(paste("religiosity", x)), data = working, Hess = TRUE)
h4.med2 <- polr(as.formula(paste("religiosity", x)), data = working2, Hess = TRUE)

## Outcomes
h4.therm <- lm(as.formula(paste("therm", x, " + religiosity")), data = working)
h4.p1 <- lm(as.formula(paste("p1.wash", x, " + religiosity")), data = working)
h4.p2 <- lm(as.formula(paste("p2.dist", x, " + religiosity")), data = working)
h4.p3 <- lm(as.formula(paste("p3.mask", x, " + religiosity")), data = working)
h4.p4 <- lm(as.formula(paste("p4.pray", x, " + religiosity")), data = working)
h4.behave <- lm(as.formula(paste("behave", x, " + religiosity")), data = working)
h4.click <- lm(as.formula(paste("clicked", x, " + religiosity")), data = working2)

## Endorsement
h4.e.therm <- mediate(h4.med, h4.therm, treat = "treat", mediator = "religiosity", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

h4.e.p1 <- mediate(h4.med, h4.p1, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

h4.e.p2 <- mediate(h4.med, h4.p2, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

h4.e.p3 <- mediate(h4.med, h4.p3, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

h4.e.p4 <- mediate(h4.med, h4.p4, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

h4.e.behave <- mediate(h4.med, h4.behave, treat = "treat", mediator = "religiosity", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

h4.e.click <- mediate(h4.med2, h4.click, treat = "treat", mediator = "religiosity", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

## Norms
h4.n.therm <- mediate(h4.med, h4.therm, treat = "treat", mediator = "religiosity", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 100)

h4.n.p1 <- mediate(h4.med, h4.p1, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Norms", robustSE = TRUE, sims = 100)

h4.n.p2 <- mediate(h4.med, h4.p2, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Norms", robustSE = TRUE, sims = 100)

h4.n.p3 <- mediate(h4.med, h4.p3, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Norms", robustSE = TRUE, sims = 100)

h4.n.p4 <- mediate(h4.med, h4.p4, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Norms", robustSE = TRUE, sims = 100)

h4.n.behave <- mediate(h4.med, h4.behave, treat = "treat", mediator = "religiosity", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Norms", robustSE = TRUE, sims = 100)

h4.n.click <- mediate(h4.med2, h4.click, treat = "treat", mediator = "religiosity", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 100)

## Mediation Plots
## Endorsement
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(h4.e.therm, main = "Thermometer")
plot(h4.e.p1, main = "Policy: Hand Washing")
plot(h4.e.p2, main = "Policy: Distancing")
plot(h4.e.p3, main = "Policy: Mask Wearing")
plot(h4.e.p4, main = "Policy: Prayer")
plot(h4.e.behave, main = "Would Click")
plot(h4.e.click, main = "Clicked")
mtext("Mediator Variable: Religiosity (Self-Reported) \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Norms
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(h4.n.therm, main = "Thermometer")
plot(h4.n.p1, main = "Policy: Hand Washing")
plot(h4.n.p2, main = "Policy: Distancing")
plot(h4.n.p3, main = "Policy: Mask Wearing")
plot(h4.n.p4, main = "Policy: Prayer")
plot(h4.n.behave, main = "Would Click")
plot(h4.n.click, main = "Clicked")
mtext("Mediator Variable: Religiosity (Self-Reported) \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)


## Religious Motivations #######################################################
rm(list=setdiff(ls(), c("working", "working2", "x")))

## Mediator
mc1 <- lm(paste("c1", x), data = working)
mc2 <- lm(paste("c2", x), data = working)
mc3 <- lm(paste("c3", x), data = working)
mc4 <- lm(paste("c4", x), data = working)
mc5 <- lm(paste("c5", x), data = working)
mc6 <- lm(paste("c6", x), data = working)

mc1.2 <- lm(paste("c1", x), data = working2)
mc2.2 <- lm(paste("c2", x), data = working2)
mc3.2 <- lm(paste("c3", x), data = working2)
mc4.2 <- lm(paste("c4", x), data = working2)
mc5.2 <- lm(paste("c5", x), data = working2)
mc6.2 <- lm(paste("c6", x), data = working2)

## Outcomes
## Motivation 1
o1c1 <- lm(paste("therm", x, " + c1"), data = working)
o2c1 <- lm(paste("p1.wash", x, " + c1"), data = working)
o3c1 <- lm(paste("p2.dist", x, " + c1"), data = working)
o4c1 <- lm(paste("p3.mask", x, " + c1"), data = working)
o5c1 <- lm(paste("p4.pray", x, " + c1"), data = working)
o6c1 <- lm(paste("behave", x, " + c1"), data = working)
o7c1 <- lm(paste("clicked", x, " + c1"), data = working2)

## Motivation 2
o1c2 <- lm(paste("therm", x, " + c2"), data = working)
o2c2 <- lm(paste("p1.wash", x, " + c2"), data = working)
o3c2 <- lm(paste("p2.dist", x, " + c2"), data = working)
o4c2 <- lm(paste("p3.mask", x, " + c2"), data = working)
o5c2 <- lm(paste("p4.pray", x, " + c2"), data = working)
o6c2 <- lm(paste("behave", x, " + c2"), data = working)
o7c2 <- lm(paste("clicked", x, " + c2"), data = working2)

## Motivation 3
o1c3 <- lm(paste("therm", x, " + c3"), data = working)
o2c3 <- lm(paste("p1.wash", x, " + c3"), data = working)
o3c3 <- lm(paste("p2.dist", x, " + c3"), data = working)
o4c3 <- lm(paste("p3.mask", x, " + c3"), data = working)
o5c3 <- lm(paste("p4.pray", x, " + c3"), data = working)
o6c3 <- lm(paste("behave", x, " + c3"), data = working)
o7c3 <- lm(paste("clicked", x, " + c3"), data = working2)

## Motivation 4
o1c4 <- lm(paste("therm", x, " + c4"), data = working)
o2c4 <- lm(paste("p1.wash", x, " + c4"), data = working)
o3c4 <- lm(paste("p2.dist", x, " + c4"), data = working)
o4c4 <- lm(paste("p3.mask", x, " + c4"), data = working)
o5c4 <- lm(paste("p4.pray", x, " + c4"), data = working)
o6c4 <- lm(paste("behave", x, " + c4"), data = working)
o7c4 <- lm(paste("clicked", x, " + c4"), data = working2)

## Motivation 5
o1c5 <- lm(paste("therm", x, " + c5"), data = working)
o2c5 <- lm(paste("p1.wash", x, " + c5"), data = working)
o3c5 <- lm(paste("p2.dist", x, " + c5"), data = working)
o4c5 <- lm(paste("p3.mask", x, " + c5"), data = working)
o5c5 <- lm(paste("p4.pray", x, " + c5"), data = working)
o6c5 <- lm(paste("behave", x, " + c5"), data = working)
o7c5 <- lm(paste("clicked", x, " + c5"), data = working2)

## Motivation 6
o1c6 <- lm(paste("therm", x, " + c6"), data = working)
o2c6 <- lm(paste("p1.wash", x, " + c6"), data = working)
o3c6 <- lm(paste("p2.dist", x, " + c6"), data = working)
o4c6 <- lm(paste("p3.mask", x, " + c6"), data = working)
o5c6 <- lm(paste("p4.pray", x, " + c6"), data = working)
o6c6 <- lm(paste("behave", x, " + c6"), data = working)
o7c6 <- lm(paste("clicked", x, " + c6"), data = working2)

## Mediation
## Motivation 1, Endorsement Condition
m1o1c1 <- mediate(mc1, o1c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o2c1 <- mediate(mc1, o2c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o3c1 <- mediate(mc1, o3c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o4c1 <- mediate(mc1, o4c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o5c1 <- mediate(mc1, o5c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o6c1 <- mediate(mc1, o6c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o7c1 <- mediate(mc1.2, o7c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

## Motivation 1, Norms Condition
m2o1c1 <- mediate(mc1, o1c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o2c1 <- mediate(mc1, o2c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o3c1 <- mediate(mc1, o3c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o4c1 <- mediate(mc1, o4c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o5c1 <- mediate(mc1, o5c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o6c1 <- mediate(mc1, o6c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o7c1 <- mediate(mc1.2, o7c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

## Motivation 2, Endorsement Condition
m1o1c2 <- mediate(mc2, o1c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o2c2 <- mediate(mc2, o2c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o3c2 <- mediate(mc2, o3c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o4c2 <- mediate(mc2, o4c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o5c2 <- mediate(mc2, o5c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o6c2 <- mediate(mc2, o6c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o7c2 <- mediate(mc2.2, o7c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

## Motivation 2, Norms Condition
m2o1c2 <- mediate(mc2, o1c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o2c2 <- mediate(mc2, o2c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o3c2 <- mediate(mc2, o3c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o4c2 <- mediate(mc2, o4c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o5c2 <- mediate(mc2, o5c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o6c2 <- mediate(mc2, o6c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o7c2 <- mediate(mc2.2, o7c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

## Motivation 3, Endorsement Condition
m1o1c3 <- mediate(mc3, o1c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o2c3 <- mediate(mc3, o2c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o3c3 <- mediate(mc3, o3c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o4c3 <- mediate(mc3, o4c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o5c3 <- mediate(mc3, o5c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o6c3 <- mediate(mc3, o6c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o7c3 <- mediate(mc3.2, o7c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

## Motivation 3, Norms Condition
m2o1c3 <- mediate(mc3, o1c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o2c3 <- mediate(mc3, o2c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o3c3 <- mediate(mc3, o3c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o4c3 <- mediate(mc3, o4c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o5c3 <- mediate(mc3, o5c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o6c3 <- mediate(mc3, o6c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o7c3 <- mediate(mc3.2, o7c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

## Motivation 4, Endorsement Condition
m1o1c4 <- mediate(mc4, o1c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o2c4 <- mediate(mc4, o2c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o3c4 <- mediate(mc4, o3c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o4c4 <- mediate(mc4, o4c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o5c4 <- mediate(mc4, o5c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o6c4 <- mediate(mc4, o6c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o674 <- mediate(mc4.2, o7c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

## Motivation 4, Norms Condition
m2o1c4 <- mediate(mc4, o1c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o2c4 <- mediate(mc4, o2c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o3c4 <- mediate(mc4, o3c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o4c4 <- mediate(mc4, o4c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o5c4 <- mediate(mc4, o5c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o6c4 <- mediate(mc4, o6c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o7c4 <- mediate(mc4.2, o7c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

## Motivation 5, Endorsement Condition
m1o1c5 <- mediate(mc5, o1c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o2c5 <- mediate(mc5, o2c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o3c5 <- mediate(mc5, o3c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o4c5 <- mediate(mc5, o4c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o5c5 <- mediate(mc5, o5c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o6c5 <- mediate(mc5, o6c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o7c5 <- mediate(mc5.2, o7c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

## Motivation 5, Norms Condition
m2o1c5 <- mediate(mc5, o1c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o2c5 <- mediate(mc5, o2c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o3c5 <- mediate(mc5, o3c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o4c5 <- mediate(mc5, o4c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o5c5 <- mediate(mc5, o5c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o6c5 <- mediate(mc5, o6c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o7c5 <- mediate(mc5.2, o7c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

## Motivation 6, Endorsement Condition
m1o1c6 <- mediate(mc6, o1c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o2c6 <- mediate(mc6, o2c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o3c6 <- mediate(mc6, o3c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o4c6 <- mediate(mc6, o4c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o5c6 <- mediate(mc6, o5c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o6c6 <- mediate(mc6, o6c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

m1o7c6 <- mediate(mc6.2, o7c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 100)

## Motivation 6, Norms Condition
m2o1c6 <- mediate(mc6, o1c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o2c6 <- mediate(mc6, o2c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o3c6 <- mediate(mc6, o3c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o4c6 <- mediate(mc6, o4c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o5c6 <- mediate(mc6, o5c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o6c6 <- mediate(mc6, o6c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

m2o7c6 <- mediate(mc6.2, o7c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 100)

## Mediation Plots
## Motivation 1, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m1o1c1, main = "Thermometer")
plot(m1o2c1, main = "Policy: Hand Washing")
plot(m1o3c1, main = "Policy: Distancing")
plot(m1o4c1, main = "Policy: Mask Wearing")
plot(m1o5c1, main = "Policy: Prayer")
plot(m1o6c1, main = "Would Click")
plot(m1o7c1, main = "Clicked")
mtext("Mediator Variable: Insurance \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 1, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m2o1c1, main = "Thermometer")
plot(m2o2c1, main = "Policy: Hand Washing")
plot(m2o3c1, main = "Policy: Distancing")
plot(m2o4c1, main = "Policy: Mask Wearing")
plot(m2o5c1, main = "Policy: Prayer")
plot(m2o6c1, main = "Would Click")
plot(m2o7c1, main = "Clicked")
mtext("Mediator Variable: Insurance \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 2, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m1o1c2, main = "Thermometer")
plot(m1o2c2, main = "Policy: Hand Washing")
plot(m1o3c2, main = "Policy: Distancing")
plot(m1o4c2, main = "Policy: Mask Wearing")
plot(m1o5c2, main = "Policy: Prayer")
plot(m1o6c2, main = "Would Click")
plot(m1o7c2, main = "Clicked")
mtext("Mediator Variable: Morality \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 2, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m2o1c2, main = "Thermometer")
plot(m2o2c2, main = "Policy: Hand Washing")
plot(m2o3c2, main = "Policy: Distancing")
plot(m2o4c2, main = "Policy: Mask Wearing")
plot(m2o5c2, main = "Policy: Prayer")
plot(m2o6c2, main = "Would Click")
plot(m2o7c2, main = "Clicked")
mtext("Mediator Variable: Morality \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 3, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m1o1c3, main = "Thermometer")
plot(m1o2c3, main = "Policy: Hand Washing")
plot(m1o3c3, main = "Policy: Distancing")
plot(m1o4c3, main = "Policy: Mask Wearing")
plot(m1o5c3, main = "Policy: Prayer")
plot(m1o6c3, main = "Would Click")
plot(m1o7c3, main = "Clicked")
mtext("Mediator Variable: Social Capital \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 3, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m2o1c3, main = "Thermometer")
plot(m2o2c3, main = "Policy: Hand Washing")
plot(m2o3c3, main = "Policy: Distancing")
plot(m2o4c3, main = "Policy: Mask Wearing")
plot(m2o5c3, main = "Policy: Prayer")
plot(m2o6c3, main = "Would Click")
plot(m2o7c3, main = "Clicked")
mtext("Mediator Variable: Social Capital \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 4, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m1o1c4, main = "Thermometer")
plot(m1o2c4, main = "Policy: Hand Washing")
plot(m1o3c4, main = "Policy: Distancing")
plot(m1o4c4, main = "Policy: Mask Wearing")
plot(m1o5c4, main = "Policy: Prayer")
plot(m1o6c4, main = "Would Click")
plot(m1o674, main = "Clicked")
mtext("Mediator Variable: Heuristics \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 4, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m2o1c4, main = "Thermometer")
plot(m2o2c4, main = "Policy: Hand Washing")
plot(m2o3c4, main = "Policy: Distancing")
plot(m2o4c4, main = "Policy: Mask Wearing")
plot(m2o5c4, main = "Policy: Prayer")
plot(m2o6c4, main = "Would Click")
plot(m2o7c4, main = "Clicked")
mtext("Mediator Variable:  Heuristics \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 5, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m1o1c5, main = "Thermometer")
plot(m1o2c5, main = "Policy: Hand Washing")
plot(m1o3c5, main = "Policy: Distancing")
plot(m1o4c5, main = "Policy: Mask Wearing")
plot(m1o5c5, main = "Policy: Prayer")
plot(m1o6c5, main = "Would Click")
plot(m1o7c5, main = "Clicked")
mtext("Mediator Variable: Tradition \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 5, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m2o1c5, main = "Thermometer")
plot(m2o2c5, main = "Policy: Hand Washing")
plot(m2o3c5, main = "Policy: Distancing")
plot(m2o4c5, main = "Policy: Mask Wearing")
plot(m2o5c5, main = "Policy: Prayer")
plot(m2o6c5, main = "Would Click")
plot(m2o7c5, main = "Clicked")
mtext("Mediator Variable: Tradition \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 6, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m1o1c6, main = "Thermometer")
plot(m1o2c6, main = "Policy: Hand Washing")
plot(m1o3c6, main = "Policy: Distancing")
plot(m1o4c6, main = "Policy: Mask Wearing")
plot(m1o5c6, main = "Policy: Prayer")
plot(m1o6c6, main = "Would Click")
plot(m1o7c6, main = "Clicked")
mtext("Mediator Variable: Social Identity \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 6, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m2o1c6, main = "Thermometer")
plot(m2o2c6, main = "Policy: Hand Washing")
plot(m2o3c6, main = "Policy: Distancing")
plot(m2o4c6, main = "Policy: Mask Wearing")
plot(m2o5c6, main = "Policy: Prayer")
plot(m2o6c6, main = "Would Click")
plot(m2o7c6, main = "Clicked")
mtext("Mediator Variable: Social Identity \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Cleaning Up
rm(list=setdiff(ls(), c("working", "x")))


## Extrinsic / Intrinsic Motivations ###########################################
## Data Cleaning
working$ei1.f <- factor(working$ei1)
working$ei2.f <- factor(working$ei2)
working$ei3.f <- factor(working$ei3)

working$ei1.f <- dplyr::recode(working$ei1.f,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Neither Agree nor Disagree",
                               "4" = "Agree",
                               "5" = "Strongly Agree") 

working$ei2.f <- dplyr::recode(working$ei2.f,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Neither Agree nor Disagree",
                               "4" = "Agree",
                               "5" = "Strongly Agree") 

working$ei3.f <- dplyr::recode(working$ei3.f,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Neither Agree nor Disagree",
                               "4" = "Agree",
                               "5" = "Strongly Agree") 

working2 <- working %>% 
  filter(behave > 4)

## Mediator
m.ei1 <- polr(as.formula(paste("ei1.f", x)), data = working, Hess = TRUE)
m.ei2 <- polr(as.formula(paste("ei2.f", x)), data = working, Hess = TRUE)
m.ei3 <- polr(as.formula(paste("ei3.f", x)), data = working, Hess = TRUE)

m2.ei1 <- polr(as.formula(paste("ei1.f", x)), data = working2, Hess = TRUE)
m2.ei2 <- polr(as.formula(paste("ei2.f", x)), data = working2, Hess = TRUE)
m2.ei3 <- polr(as.formula(paste("ei3.f", x)), data = working2, Hess = TRUE)

## Outcomes
## EI 1 - Religion Offers Comfort
ei1.therm <- lm(as.formula(paste("therm", x, " + ei1.f")), data = working)
ei1.p1 <- lm(as.formula(paste("p1.wash", x, " + ei1.f")), data = working)
ei1.p2 <- lm(as.formula(paste("p2.dist", x, " + ei1.f")), data = working)
ei1.p3 <- lm(as.formula(paste("p3.mask", x, " + ei1.f")), data = working)
ei1.p4 <- lm(as.formula(paste("p4.pray", x, " + ei1.f")), data = working)
ei1.behave <- lm(as.formula(paste("behave", x, " + ei1.f")), data = working)
ei1.click <- lm(as.formula(paste("clicked", x, " + ei1.f")), data = working2)

## EI 2 - Religious Approach to Life
ei2.therm <- lm(as.formula(paste("therm", x, " + ei2.f")), data = working)
ei2.p1 <- lm(as.formula(paste("p1.wash", x, " + ei2.f")), data = working)
ei2.p2 <- lm(as.formula(paste("p2.dist", x, " + ei2.f")), data = working)
ei2.p3 <- lm(as.formula(paste("p3.mask", x, " + ei2.f")), data = working)
ei2.p4 <- lm(as.formula(paste("p4.pray", x, " + ei2.f")), data = working)
ei2.behave <- lm(as.formula(paste("behave", x, " + ei2.f")), data = working)
ei2.click <- lm(as.formula(paste("clicked", x, " + ei2.f")), data = working2)

## EI 3 - Enjoys Social Aspects of Church
ei3.therm <- lm(as.formula(paste("therm", x, " + ei3.f")), data = working)
ei3.p1 <- lm(as.formula(paste("p1.wash", x, " + ei3.f")), data = working)
ei3.p2 <- lm(as.formula(paste("p2.dist", x, " + ei3.f")), data = working)
ei3.p3 <- lm(as.formula(paste("p3.mask", x, " + ei3.f")), data = working)
ei3.p4 <- lm(as.formula(paste("p4.pray", x, " + ei3.f")), data = working)
ei3.behave <- lm(as.formula(paste("behave", x, " + ei3.f")), data = working)
ei3.click <- lm(as.formula(paste("clicked", x, " + ei3.f")), data = working2)

## Mediation
## EI 1,  Endorsement Condition
ei1.e.therm <- mediate(m.ei1, ei1.therm, treat = "treat", mediator = "ei1.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei1.e.p1 <- mediate(m.ei1, ei1.p1, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei1.e.p2 <- mediate(m.ei1, ei1.p2, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei1.e.p3 <- mediate(m.ei1, ei1.p3, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei1.e.p4 <- mediate(m.ei1, ei1.p4, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei1.e.behave <- mediate(m.ei1, ei1.behave, treat = "treat", mediator = "ei1.f", 
                        dropobs = TRUE, control.value = "Control", 
                        treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei1.e.click <- mediate(m2.ei1, ei1.click, treat = "treat", mediator = "ei1.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

## EI 1, Norms Condition
ei1.n.therm <- mediate(m.ei1, ei1.therm, treat = "treat", mediator = "ei1.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei1.n.p1 <- mediate(m.ei1, ei1.p1, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei1.n.p2 <- mediate(m.ei1, ei1.p2, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei1.n.p3 <- mediate(m.ei1, ei1.p3, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei1.n.p4 <- mediate(m.ei1, ei1.p4, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei1.n.behave <- mediate(m.ei1, ei1.behave, treat = "treat", mediator = "ei1.f", 
                        dropobs = TRUE, control.value = "Control", 
                        treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei1.n.click <- mediate(m2.ei1, ei1.click, treat = "treat", mediator = "ei1.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Norms", robustSE = TRUE, sims = 100)

## EI2, Endorsement Condition
ei2.e.therm <- mediate(m.ei2, ei2.therm, treat = "treat", mediator = "ei2.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei2.e.p1 <- mediate(m.ei2, ei2.p1, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei2.e.p2 <- mediate(m.ei2, ei2.p2, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei2.e.p3 <- mediate(m.ei2, ei2.p3, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei2.e.p4 <- mediate(m.ei2, ei2.p4, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei2.e.behave <- mediate(m.ei2, ei2.behave, treat = "treat", mediator = "ei2.f", 
                        dropobs = TRUE, control.value = "Control", 
                        treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei2.e.click <- mediate(m2.ei2, ei2.click, treat = "treat", mediator = "ei2.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

## EI 2, Norms Condition
ei2.n.therm <- mediate(m.ei2, ei2.therm, treat = "treat", mediator = "ei2.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei2.n.p1 <- mediate(m.ei2, ei2.p1, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei2.n.p2 <- mediate(m.ei2, ei2.p2, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei2.n.p3 <- mediate(m.ei2, ei2.p3, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei2.n.p4 <- mediate(m.ei2, ei2.p4, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei2.n.behave <- mediate(m.ei2, ei2.behave, treat = "treat", mediator = "ei2.f", 
                        dropobs = TRUE, control.value = "Control", 
                        treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei2.n.click <- mediate(m2.ei2, ei2.click, treat = "treat", mediator = "ei2.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Norms", robustSE = TRUE, sims = 100)

## EI 3, Endorsement Condition
ei3.e.therm <- mediate(m.ei3, ei3.therm, treat = "treat", mediator = "ei3.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei3.e.p1 <- mediate(m.ei3, ei3.p1, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei3.e.p2 <- mediate(m.ei3, ei3.p2, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei3.e.p3 <- mediate(m.ei3, ei3.p3, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei3.e.p4 <- mediate(m.ei3, ei3.p4, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei3.e.behave <- mediate(m.ei3, ei3.behave, treat = "treat", mediator = "ei3.f", 
                        dropobs = TRUE, control.value = "Control", 
                        treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

ei3.e.click <- mediate(m2.ei3, ei3.click, treat = "treat", mediator = "ei3.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

## EI 3, Norms Condition
ei3.n.therm <- mediate(m.ei3, ei3.therm, treat = "treat", mediator = "ei3.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei3.n.p1 <- mediate(m.ei3, ei3.p1, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei3.n.p2 <- mediate(m.ei3, ei3.p2, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei3.n.p3 <- mediate(m.ei3, ei3.p3, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei3.n.p4 <- mediate(m.ei3, ei3.p4, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei3.n.behave <- mediate(m.ei3, ei3.behave, treat = "treat", mediator = "ei3.f", 
                        dropobs = TRUE, control.value = "Control", 
                        treat.value  = "Norms", robustSE = TRUE, sims = 100)

ei3.n.click <- mediate(m2.ei3, ei3.click, treat = "treat", mediator = "ei3.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Norms", robustSE = TRUE, sims = 100)

## Mediation Plots
## EI 1, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(ei1.e.therm, main = "Thermometer")
plot(ei1.e.p1, main = "Policy: Hand Washing")
plot(ei1.e.p2, main = "Policy: Distancing")
plot(ei1.e.p3, main = "Policy: Mask Wearing")
plot(ei1.e.p4, main = "Policy: Prayer")
plot(ei1.e.behave, main = "Would Click")
plot(ei1.e.click, main = "Clicked")
mtext("Mediator Variable: Religion Offers Comfort \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## EI 1, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(ei1.n.therm, main = "Thermometer")
plot(ei1.n.p1, main = "Policy: Hand Washing")
plot(ei1.n.p2, main = "Policy: Distancing")
plot(ei1.n.p3, main = "Policy: Mask Wearing")
plot(ei1.n.p4, main = "Policy: Prayer")
plot(ei1.n.behave, main = "Would Click")
plot(ei1.n.click, main = "Clicked")
mtext("Mediator Variable: Religion Offers Comfort \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## EI 2, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(ei2.e.therm, main = "Thermometer")
plot(ei2.e.p1, main = "Policy: Hand Washing")
plot(ei2.e.p2, main = "Policy: Distancing")
plot(ei2.e.p3, main = "Policy: Mask Wearing")
plot(ei2.e.p4, main = "Policy: Prayer")
plot(ei2.e.behave, main = "Would Click")
plot(ei2.e.click, main = "Clicked")
mtext("Mediator Variable: Religious Approach to Life \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## EI 2, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(ei2.n.therm, main = "Thermometer")
plot(ei2.n.p1, main = "Policy: Hand Washing")
plot(ei2.n.p2, main = "Policy: Distancing")
plot(ei2.n.p3, main = "Policy: Mask Wearing")
plot(ei2.n.p4, main = "Policy: Prayer")
plot(ei2.n.behave, main = "Would Click")
plot(ei2.n.click, main = "Clicked")
mtext("Mediator Variable: Religious Approach to Life \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## EI 3, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(ei3.e.therm, main = "Thermometer")
plot(ei3.e.p1, main = "Policy: Hand Washing")
plot(ei3.e.p2, main = "Policy: Distancing")
plot(ei3.e.p3, main = "Policy: Mask Wearing")
plot(ei3.e.p4, main = "Policy: Prayer")
plot(ei3.e.behave, main = "Would Click")
plot(ei3.e.click, main = "Clicked")
mtext("Mediator Variable: Enjoys Social Aspects of Church \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## EI 3, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(ei3.n.therm, main = "Thermometer")
plot(ei3.n.p1, main = "Policy: Hand Washing")
plot(ei3.n.p2, main = "Policy: Distancing")
plot(ei3.n.p3, main = "Policy: Mask Wearing")
plot(ei3.n.p4, main = "Policy: Prayer")
plot(ei3.n.behave, main = "Would Click")
plot(ei3.n.click, main = "Clicked")
mtext("Mediator Variable: Enjoys Social Aspects of Church \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Cleaning Up
rm(list=setdiff(ls(), c("working", "x")))


## Religiosity Index ###########################################################
working$rel1x <- factor(working$rel1x)
working$rel2x <- factor(working$rel2x)
working$rel3x <- factor(working$rel3x)

working$rel1x <- dplyr::recode(working$rel1x,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Somewhat Disagree",
                               "4" = "Somewhat Agree",
                               "5" = "Agree",
                               "6" = "Strongly Agree")

working$rel2x <- dplyr::recode(working$rel2x,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Somewhat Disagree",
                               "4" = "Somewhat Agree",
                               "5" = "Agree",
                               "6" = "Strongly Agree")

working$rel3x <- dplyr::recode(working$rel3x,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Somewhat Disagree",
                               "4" = "Somewhat Agree",
                               "5" = "Agree",
                               "6" = "Strongly Agree")

working2 <- working %>% 
  filter(behave > 4)

## Mediator
m.rel1x <- polr(as.formula(paste("rel1x", x)), data = working, Hess = TRUE)
m.rel2x <- polr(as.formula(paste("rel2x", x)), data = working, Hess = TRUE)
m.rel3x <- polr(as.formula(paste("rel3x", x)), data = working, Hess = TRUE)

m2.rel1x <- polr(as.formula(paste("rel1x", x)), data = working2, Hess = TRUE)
m2.rel2x <- polr(as.formula(paste("rel2x", x)), data = working2, Hess = TRUE)
m2.rel3x <- polr(as.formula(paste("rel3x", x)), data = working2, Hess = TRUE)

## Outcomes
## Religiosity Index 1 - Regularly Attends Church
rel1x.therm <- lm(as.formula(paste("therm", x, " + rel1x")), data = working)
rel1x.p1 <- lm(as.formula(paste("p1.wash", x, " + rel1x")), data = working)
rel1x.p2 <- lm(as.formula(paste("p2.dist", x, " + rel1x")), data = working)
rel1x.p3 <- lm(as.formula(paste("p3.mask", x, " + rel1x")), data = working)
rel1x.p4 <- lm(as.formula(paste("p4.pray", x, " + rel1x")), data = working)
rel1x.behave <- lm(as.formula(paste("behave", x, " + rel1x")), data = working)
rel1x.click <- lm(as.formula(paste("clicked", x, " + rel1x")), data = working2)

## Religiosity Index 2 - Spiritual Values
rel2x.therm <- lm(as.formula(paste("therm", x, " + rel2x")), data = working)
rel2x.p1 <- lm(as.formula(paste("p1.wash", x, " + rel2x")), data = working)
rel2x.p2 <- lm(as.formula(paste("p2.dist", x, " + rel2x")), data = working)
rel2x.p3 <- lm(as.formula(paste("p3.mask", x, " + rel2x")), data = working)
rel2x.p4 <- lm(as.formula(paste("p4.pray", x, " + rel2x")), data = working)
rel2x.behave <- lm(as.formula(paste("behave", x, " + rel2x")), data = working)
rel2x.click <- lm(as.formula(paste("clicked", x, " + rel2x")), data = working2)

## Religiosity Index 3 - Enjoys Social Aspects
rel3x.therm <- lm(as.formula(paste("therm", x, " + rel3x")), data = working)
rel3x.p1 <- lm(as.formula(paste("p1.wash", x, " + rel3x")), data = working)
rel3x.p2 <- lm(as.formula(paste("p2.dist", x, " + rel3x")), data = working)
rel3x.p3 <- lm(as.formula(paste("p3.mask", x, " + rel3x")), data = working)
rel3x.p4 <- lm(as.formula(paste("p4.pray", x, " + rel3x")), data = working)
rel3x.behave <- lm(as.formula(paste("behave", x, " + rel3x")), data = working)
rel3x.click <- lm(as.formula(paste("clicked", x, " + rel3x")), data = working2)

## Mediation
## Rel Index 1, Endorsement Condition
rel1x.e.therm <- mediate(m.rel1x, rel1x.therm, treat = "treat", mediator = "rel1x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel1x.e.p1 <- mediate(m.rel1x, rel1x.p1, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel1x.e.p2 <- mediate(m.rel1x, rel1x.p2, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel1x.e.p3 <- mediate(m.rel1x, rel1x.p3, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel1x.e.p4 <- mediate(m.rel1x, rel1x.p4, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel1x.e.behave <- mediate(m.rel1x, rel1x.behave, treat = "treat", mediator = "rel1x", 
                          dropobs = TRUE, control.value = "Control", 
                          treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel1x.e.click <- mediate(m2.rel1x, rel1x.click, treat = "treat", mediator = "rel1x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

## Rel Index 1, Norms Condition
rel1x.n.therm <- mediate(m.rel1x, rel1x.therm, treat = "treat", mediator = "rel1x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel1x.n.p1 <- mediate(m.rel1x, rel1x.p1, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel1x.n.p2 <- mediate(m.rel1x, rel1x.p2, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel1x.n.p3 <- mediate(m.rel1x, rel1x.p3, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel1x.n.p4 <- mediate(m.rel1x, rel1x.p4, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel1x.n.behave <- mediate(m.rel1x, rel1x.behave, treat = "treat", mediator = "rel1x", 
                          dropobs = TRUE, control.value = "Control", 
                          treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel1x.n.click <- mediate(m2.rel1x, rel1x.click, treat = "treat", mediator = "rel1x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Norms", robustSE = TRUE, sims = 100)

## Rel Index 2, Endorsement Condition
rel2x.e.therm <- mediate(m.rel2x, rel2x.therm, treat = "treat", mediator = "rel2x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel2x.e.p1 <- mediate(m.rel2x, rel2x.p1, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel2x.e.p2 <- mediate(m.rel2x, rel2x.p2, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel2x.e.p3 <- mediate(m.rel2x, rel2x.p3, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel2x.e.p4 <- mediate(m.rel2x, rel2x.p4, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel2x.e.behave <- mediate(m.rel2x, rel2x.behave, treat = "treat", mediator = "rel2x", 
                          dropobs = TRUE, control.value = "Control", 
                          treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel2x.e.click <- mediate(m2.rel2x, rel2x.click, treat = "treat", mediator = "rel2x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

## Rel Index 2, Norms Condition
rel2x.n.therm <- mediate(m.rel2x, rel2x.therm, treat = "treat", mediator = "rel2x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel2x.n.p1 <- mediate(m.rel2x, rel2x.p1, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel2x.n.p2 <- mediate(m.rel2x, rel2x.p2, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel2x.n.p3 <- mediate(m.rel2x, rel2x.p3, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel2x.n.p4 <- mediate(m.rel2x, rel2x.p4, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel2x.n.behave <- mediate(m.rel2x, rel2x.behave, treat = "treat", mediator = "rel2x", 
                          dropobs = TRUE, control.value = "Control", 
                          treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel2x.n.click <- mediate(m2.rel2x, rel2x.click, treat = "treat", mediator = "rel2x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Norms", robustSE = TRUE, sims = 100)

## Rel Index 3, Endorsement Condition
rel3x.e.therm <- mediate(m.rel3x, rel3x.therm, treat = "treat", mediator = "rel3x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel3x.e.p1 <- mediate(m.rel3x, rel3x.p1, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel3x.e.p2 <- mediate(m.rel3x, rel3x.p2, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel3x.e.p3 <- mediate(m.rel3x, rel3x.p3, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel3x.e.p4 <- mediate(m.rel3x, rel3x.p4, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel3x.e.behave <- mediate(m.rel3x, rel3x.behave, treat = "treat", mediator = "rel3x", 
                          dropobs = TRUE, control.value = "Control", 
                          treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

rel3x.e.click <- mediate(m2.rel3x, rel3x.click, treat = "treat", mediator = "rel3x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Endorsement", robustSE = TRUE, sims = 100)

## Rel Index 3, Norms Condition
rel3x.n.therm <- mediate(m.rel3x, rel3x.therm, treat = "treat", mediator = "rel3x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel3x.n.p1 <- mediate(m.rel3x, rel3x.p1, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel3x.n.p2 <- mediate(m.rel3x, rel3x.p2, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel3x.n.p3 <- mediate(m.rel3x, rel3x.p3, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel3x.n.p4 <- mediate(m.rel3x, rel3x.p4, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel3x.n.behave <- mediate(m.rel3x, rel3x.behave, treat = "treat", mediator = "rel3x", 
                          dropobs = TRUE, control.value = "Control", 
                          treat.value  = "Norms", robustSE = TRUE, sims = 100)

rel3x.n.click <- mediate(m2.rel3x, rel3x.click, treat = "treat", mediator = "rel3x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Norms", robustSE = TRUE, sims = 100)

## Mediation Plots (rel1x)
## Rel Index 1, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(rel1x.e.therm, main = "Thermometer")
plot(rel1x.e.p1, main = "Policy: Hand Washing")
plot(rel1x.e.p2, main = "Policy: Distancing")
plot(rel1x.e.p3, main = "Policy: Mask Wearing")
plot(rel1x.e.p4, main = "Policy: Prayer")
plot(rel1x.e.behave, main = "Would Click")
plot(rel1x.e.click, main = "Clicked")
mtext("Mediator Variable: Regularly Attends Church \n Endorsement Condition", outer=TRUE,  side = 3, cex=1, line=-0.5)

## Rel Index 1, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(rel1x.n.therm, main = "Thermometer")
plot(rel1x.n.p1, main = "Policy: Hand Washing")
plot(rel1x.n.p2, main = "Policy: Distancing")
plot(rel1x.n.p3, main = "Policy: Mask Wearing")
plot(rel1x.n.p4, main = "Policy: Prayer")
plot(rel1x.n.behave, main = "Would Click")
plot(rel1x.n.click, main = "Clicked")
mtext("Mediator Variable: Regularly Attends Church \n Norms Condition", outer=TRUE,  side = 3, cex=1, line=-0.5)

## Rel Index 2, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(rel2x.e.therm, main = "Thermometer")
plot(rel2x.e.p1, main = "Policy: Hand Washing")
plot(rel2x.e.p2, main = "Policy: Distancing")
plot(rel2x.e.p3, main = "Policy: Mask Wearing")
plot(rel2x.e.p4, main = "Policy: Prayer")
plot(rel2x.e.behave, main = "Would Click")
plot(rel2x.e.click, main = "Clicked")
mtext("Mediator Variable: Spiritual Values \n Endorsement Condition", outer=TRUE,  side = 3, cex=1, line=-0.5)

## Rel Index 2, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(rel2x.n.therm, main = "Thermometer")
plot(rel2x.n.p1, main = "Policy: Hand Washing")
plot(rel2x.n.p2, main = "Policy: Distancing")
plot(rel2x.n.p3, main = "Policy: Mask Wearing")
plot(rel2x.n.p4, main = "Policy: Prayer")
plot(rel2x.n.behave, main = "Would Click")
plot(rel2x.n.click, main = "Clicked")
mtext("Mediator Variable: Spiritual Values \n Norms Condition", outer=TRUE,  side = 3, cex=1, line=-0.5)

## Rel Index 3, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(rel3x.e.therm, main = "Thermometer")
plot(rel3x.e.p1, main = "Policy: Hand Washing")
plot(rel3x.e.p2, main = "Policy: Distancing")
plot(rel3x.e.p3, main = "Policy: Mask Wearing")
plot(rel3x.e.p4, main = "Policy: Prayer")
plot(rel3x.e.behave, main = "Would Click")
plot(rel3x.e.click, main = "Clicked")
mtext("Mediator Variable: Enjoys Social Aspects of Church \n Endorsement Condition", outer=TRUE,  side = 3, cex=1, line=-0.5)

## Rel Index 3, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(rel3x.n.therm, main = "Thermometer")
plot(rel3x.n.p1, main = "Policy: Hand Washing")
plot(rel3x.n.p2, main = "Policy: Distancing")
plot(rel3x.n.p3, main = "Policy: Mask Wearing")
plot(rel3x.n.p4, main = "Policy: Prayer")
plot(rel3x.n.behave, main = "Would Click")
plot(rel3x.n.click, main = "Clicked")
mtext("Mediator Variable: Enjoys Social Aspects of Church \n Norms Condition", outer=TRUE,  side = 3, cex=1, line=-0.5)

## Cleaning Up
rm(list=setdiff(ls(), c("working", "working2")))


# Interaction Models ###########################################################
## Religiosity Mechanism #######################################################
## Data Cleaning
x <- " ~ treat:religiosity + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"

## Interaction Models
int.therm <- lm(as.formula(paste("therm", x)), data = working)
int.p1 <- lm(as.formula(paste("p1.wash", x)), data = working)
int.p2 <- lm(as.formula(paste("p2.dist", x)), data = working)
int.p3 <- lm(as.formula(paste("p3.mask", x)), data = working)
int.p4 <- lm(as.formula(paste("p4.pray", x)), data = working)
int.behave <- lm(as.formula(paste("behave", x)), data = working)
int.click <- lm(as.formula(paste("clicked", x)), data = working2)

## OLS Tables
int.mod <- list("Thermometer" = int.therm, 
                "Policy: Hand Washing" = int.p1, 
                "Policy: Distancing" = int.p2, 
                "Policy: Mask Wearing" = int.p3, 
                "Policy: Prayer" = int.p4, 
                "Would Click" = int.behave, 
                "Clicked" = int.click)

cm <- c("(Intercept)" =	"Intercept",
        "treatNorms:religiosityAnti-Religious" =	"Norms : Anti-Religious",
        "treatNorms:religiosityNot at all Religious" =	"Norms : Not at all Religious",
        "treatNorms:religiositySlightly Religious" =	"Norms : Slightly Religious",
        "treatNorms:religiosityModerately Religious" =	"Norms : Moderately Religious",
        "treatNorms:religiosityVery Religious" =	"Norms : Very Religious",        
        "treatEndorsement:religiosityAnti-Religious" =	"Endorsement : Anti-Religious",
        "treatEndorsement:religiosityNot at all Religious" =	"Endorsement : Not at all Religious",
        "treatEndorsement:religiositySlightly Religious" =	"Endorsement : Slightly Religious",
        "treatEndorsement:religiosityModerately Religious" =	"Endorsement : Moderately Religious",
        "treatEndorsement:religiosityVery Religious" =	"Endorsement : Very Religious",
        "treatControl:religiosityAnti-Religious" =	"Control : Anti-Religious",
        "treatControl:religiosityNot at all Religious" =	"Control : Not at all Religious",        
        "treatControl:religiositySlightly Religious" =	"Control : Slightly Religious",
        "treatControl:religiosityModerately Religious" =	"Control : Moderately Religious",
        "treatControl:religiosityVery Religious" =	"Control : Very Religious")

modelsummary(int.mod,
             output = "latex",
             title = "Interaction Model: Treatment x Religiosity (Self-Reported)",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Margins Plot
p.int.therm <- ggpredict(int.therm, terms = c("treat", "religiosity"))
p.int.p1 <- ggpredict(int.p1, terms = c("treat", "religiosity"))
p.int.p2 <- ggpredict(int.p2, terms = c("treat", "religiosity"))
p.int.p3 <- ggpredict(int.p3, terms = c("treat", "religiosity"))
p.int.p4 <- ggpredict(int.p4, terms = c("treat", "religiosity"))
p.int.behave <- ggpredict(int.behave, terms = c("treat", "religiosity"))
p.int.click <- ggpredict(int.click, terms = c("treat", "religiosity"))

p.int.therm$model <- "Thermometer"
p.int.p1$model <- "Policy: Hand Washing"
p.int.p2$model <- "Policy: Distancing"
p.int.p3$model <- "Policy: Mask Wearing"
p.int.p4$model <- "Policy: Prayer"
p.int.behave$model <- "Would Click"
p.int.click$model <- "Clicked"

p.int.all <- rbind(p.int.therm, p.int.p1, p.int.p2, p.int.p3, p.int.p4, p.int.behave, p.int.click)

p.int.all$model <- factor(p.int.all$model, 
                          levels = c("Thermometer", 
                                     "Policy: Hand Washing", 
                                     "Policy: Distancing", 
                                     "Policy: Mask Wearing", 
                                     "Policy: Prayer", 
                                     "Would Click", 
                                     "Clicked"), 
                          ordered = TRUE)

p.int.all$x <- factor(p.int.all$x, 
                      levels = c("Norms", 
                                 "Endorsement", 
                                 "Control") ,
                      ordered = TRUE)

ggplot(p.int.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Religiosity (Self-Reported)",
       x  = "Condition", y = "Predicted Value")

## Cleaning Up
rm(list=setdiff(ls(), c("working")))

## Religious Motivations #######################################################
## Create a factor variable for the 6 religious motivations
as.data.frame(colnames(working))

working2 <- working
working2$m6 <- factor(apply(working2[,51:56], 1, function(x) names(x)[x == 1]), 
                      levels = names(working2[,51:56]))

working3 <- working2 %>%
  filter(behave > 4)

working2$treat <- factor(working2$treat, levels = c("Norms", "Endorsement", "Control"))
working3$treat <- factor(working3$treat, levels = c("Norms", "Endorsement", "Control"))

## Note: Tables in previous version of SI present the same results but changed...
## ...the omitted (reference) category
## To change the reference category for the interaction term:
# working2$treat <- relevel(working2$treat, ref = "Control") # Change Category Here
# working3$treat <- relevel(working3$treat, ref = "Control") # Change Category Here
# working2$m6 <- relevel(working2$m6, ref = "c1") # Change Category Here
# working3$m6 <- relevel(working3$m6, ref = "c1") # Change Category Here

x <- "~ m6:treat + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"

## Interaction Models
int.therm <- lm(as.formula(paste("therm", x)), data = working2)
int.p1 <- lm(as.formula(paste("p1.wash", x)), data = working2)
int.p2 <- lm(as.formula(paste("p2.dist", x)), data = working2)
int.p3 <- lm(as.formula(paste("p3.mask", x)), data = working2)
int.p4 <- lm(as.formula(paste("p4.pray", x)), data = working2)
int.behave <- lm(as.formula(paste("behave", x)), data = working2)
int.click <- lm(as.formula(paste("clicked", x)), data = working3)

## OLS Tables
m6.mod <- list("Thermometer" = int.therm, 
               "Policy: Hand Washing" = int.p1, 
               "Policy: Distancing" = int.p2, 
               "Policy: Mask Wearing" = int.p3, 
               "Policy: Prayer" = int.p4, 
               "Would Click" = int.behave, 
               "Clicked" = int.click)

cm <- c("(Intercept)" = "Intercept",
        "m6c1:treatNorms" = "Norms : Insurance",
        "m6c2:treatNorms" = "Norms : Morality",
        "m6c3:treatNorms" = "Norms : Social Capital",
        "m6c4:treatNorms" = "Norms : Heuristics",
        "m6c5:treatNorms" = "Norms : Tradition",
        "m6c6:treatNorms" = "Norms : Social Identity",
        "m6c1:treatEndorsement" = "Endorsement : Insurance",
        "m6c2:treatEndorsement" = "Endorsement : Morality",
        "m6c3:treatEndorsement" = "Endorsement : Social Capital",
        "m6c4:treatEndorsement" = "Endorsement : Heuristics",
        "m6c5:treatEndorsement" = "Endorsement : Tradition",
        "m6c6:treatEndorsement" = "Endorsement : Social Identity",
        "m6c1:treatControl" = "Control : Insurance",
        "m6c2:treatControl" = "Control : Morality",
        "m6c3:treatControl" = "Control : Social Capital",
        "m6c4:treatControl" = "Control : Heuristics",
        "m6c5:treatControl" = "Control : Tradition",
        "m6c6:treatControl" = "Control : Social Identity")

modelsummary(m6.mod,
             output = "latex",
             title = "Interaction Model: Treatment x Religious Motivations, 
             Omitted Category = Control x Social Identity",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Margins Plot
m.therm <- ggpredict(int.therm, terms = c("treat", "m6"))
m.therm$model <- "Thermometer"

m.p1 <- ggpredict(int.p1, terms = c("treat", "m6"))
m.p1$model <- "Policy: Hand Washing"

m.p2 <- ggpredict(int.p2, terms = c("treat", "m6"))
m.p2$model <- "Policy: Distancing"

m.p3 <- ggpredict(int.p3, terms = c("treat", "m6"))
m.p3$model <- "Policy: Mask Wearing"

m.p4 <- ggpredict(int.p4, terms = c("treat", "m6"))
m.p4$model <- "Policy: Prayer"

m.behave <- ggpredict(int.behave, terms = c("treat", "m6"))
m.behave$model <- "Would Click"

m.click <- ggpredict(int.click, terms = c("treat", "m6"))
m.click$model <- "Clicked"

m.combined <- rbind(m.therm, m.p1, m.p2, m.p3, m.p4, m.behave, m.click)

m.combined$model <- factor(m.combined$model, 
                           levels = c("Thermometer", 
                                      "Policy: Hand Washing", 
                                      "Policy: Distancing", 
                                      "Policy: Mask Wearing", 
                                      "Policy: Prayer", 
                                      "Would Click", 
                                      "Clicked"), 
                           ordered = TRUE)

m.combined$group <- dplyr::recode(m.combined$group,
                                  "c1" = "Insurance",
                                  "c2" = "Morality",
                                  "c3" = "Social Capital",
                                  "c4" = "Heuristics",
                                  "c5" = "Tradition",
                                  "c6" = "Social Identity")

ggplot(m.combined, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.5, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Religious Motivations",
       x  = "Condition", y = "Predicted Value")

## Cleaning Up
rm(list=setdiff(ls(), c("working", "working2")))

## Extrinsic / Intrinsic Motivation ############################################
working2 <- working %>% 
  filter(behave > 4)

iv1 <- " ~ treat:ei1.f + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"
iv2 <- " ~ treat:ei2.f + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"
iv3 <- " ~ treat:ei3.f + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"

## Interaction Models & Tables - EI 1 - Religion Offers Comfort
i1.therm <- lm(as.formula(paste("therm", iv1)), data = working)
i1.p1 <- lm(as.formula(paste("p1.wash", iv1)), data = working)
i1.p2 <- lm(as.formula(paste("p2.dist", iv1)), data = working)
i1.p3 <- lm(as.formula(paste("p3.mask", iv1)), data = working)
i1.p4 <- lm(as.formula(paste("p4.pray", iv1)), data = working)
i1.behave <- lm(as.formula(paste("behave", iv1)), data = working)
i1.click <- lm(as.formula(paste("clicked", iv1)), data = working2)

i1 <- list("Thermometer" = i1.therm, 
           "Policy: Hand Washing" = i1.p1, 
           "Policy: Distancing" = i1.p2, 
           "Policy: Mask Wearing" = i1.p3, 
           "Policy: Prayer" = i1.p4, 
           "Would Click" = i1.behave, 
           "Clicked" = i1.click)

cm <- c("(Intercept)" = "Intercept",
        "treatNorms:ei1.fStrongly Disagree" = "Norms : Strongly Disagree",
        "treatNorms:ei1.fDisagree" = "Norms : Disagree",
        "treatNorms:ei1.fNeither Agree nor Disagree" = "Norms : Neither Agree nor Disagree",
        "treatNorms:ei1.fAgree" = "Norms : Agree",
        "treatNorms:ei1.fStrongly Agree" = "Norms : Stongly Agree",
        "treatEndorsement:ei1.fStrongly Disagree" = "Endorsement : Strongly Disagree",
        "treatEndorsement:ei1.fDisagree" = "Endorsement : Disagree",
        "treatEndorsement:ei1.fNeither Agree nor Disagree" = "Endorsement : Neither Agree nor Disagree",
        "treatEndorsement:ei1.fAgree" = "Endorsement : Agree",
        "treatEndorsement:ei1.fStrongly Agree" = "Endorsement : Stongly Agree",
        "treatControl:ei1.fStrongly Disagree" = "Control : Strongly Disagree",
        "treatControl:ei1.fDisagree" = "Control : Disagree",        
        "treatControl:ei1.fNeither Agree nor Disagree" = "Control : Neither Agree nor Disagree",
        "treatControl:ei1.fAgree" = "Control : Agree",
        "treatControl:ei1.fStrongly Agree" = "Control : Stongly Agree")

modelsummary(i1,
             output = "latex",
             title = "Interaction Model: Treatment x Extrinsic / Intrinsic Motivations — Religion Offers Comfort",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Interaction Models & Tables - EI 2 - Religious Approach to Life
i2.therm <- lm(as.formula(paste("therm", iv2)), data = working)
i2.p1 <- lm(as.formula(paste("p1.wash", iv2)), data = working)
i2.p2 <- lm(as.formula(paste("p2.dist", iv2)), data = working)
i2.p3 <- lm(as.formula(paste("p3.mask", iv2)), data = working)
i2.p4 <- lm(as.formula(paste("p4.pray", iv2)), data = working)
i2.behave <- lm(as.formula(paste("behave", iv2)), data = working)
i2.click <- lm(as.formula(paste("clicked", iv2)), data = working2)

i2 <- list("Thermometer" = i2.therm, 
           "Policy: Hand Washing" = i2.p1, 
           "Policy: Distancing" = i2.p2, 
           "Policy: Mask Wearing" = i2.p3, 
           "Policy: Prayer" = i2.p4, 
           "Would Click" = i2.behave, 
           "Clicked" = i2.click)

cm <- c("(Intercept)" = "Intercept",
        "treatNorms:ei2.fStrongly Disagree" = "Norms : Strongly Disagree",
        "treatNorms:ei2.fDisagree" = "Norms : Disagree",
        "treatNorms:ei2.fNeither Agree nor Disagree" = "Norms : Neither Agree nor Disagree",
        "treatNorms:ei2.fAgree" = "Norms : Agree",
        "treatNorms:ei2.fStrongly Agree" = "Norms : Stongly Agree",
        "treatEndorsement:ei2.fStrongly Disagree" = "Endorsement : Strongly Disagree",
        "treatEndorsement:ei2.fDisagree" = "Endorsement : Disagree",
        "treatEndorsement:ei2.fNeither Agree nor Disagree" = "Endorsement : Neither Agree nor Disagree",
        "treatEndorsement:ei2.fAgree" = "Endorsement : Agree",
        "treatEndorsement:ei2.fStrongly Agree" = "Endorsement : Stongly Agree",
        "treatControl:ei2.fStrongly Disagree" = "Control : Strongly Disagree",
        "treatControl:ei2.fDisagree" = "Control : Disagree",        
        "treatControl:ei2.fNeither Agree nor Disagree" = "Control : Neither Agree nor Disagree",
        "treatControl:ei2.fAgree" = "Control : Agree",
        "treatControl:ei2.fStrongly Agree" = "Control : Stongly Agree") 

modelsummary(i2,
             output = "latex",
             title = "Interaction Model: Treatment x Extrinsic / Intrinsic Motivations — Takes Religious Approach to Life",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Interaction Models & Tables - EI 3 - Enjoys Social Aspects of Church
i3.therm <- lm(as.formula(paste("therm", iv3)), data = working)
i3.p1 <- lm(as.formula(paste("p1.wash", iv3)), data = working)
i3.p2 <- lm(as.formula(paste("p2.dist", iv3)), data = working)
i3.p3 <- lm(as.formula(paste("p3.mask", iv3)), data = working)
i3.p4 <- lm(as.formula(paste("p4.pray", iv3)), data = working)
i3.behave <- lm(as.formula(paste("behave", iv3)), data = working)
i3.click <- lm(as.formula(paste("clicked", iv3)), data = working2)

i3 <- list("Thermometer" = i3.therm, 
           "Policy: Hand Washing" = i3.p1, 
           "Policy: Distancing" = i3.p2, 
           "Policy: Mask Wearing" = i3.p3, 
           "Policy: Prayer" = i3.p4, 
           "Would Click" = i3.behave, 
           "Clicked" = i3.click)

cm <- c("(Intercept)" = "Intercept",
        "treatNorms:ei3.fStrongly Disagree" = "Norms : Strongly Disagree",
        "treatNorms:ei3.fDisagree" = "Norms : Disagree",
        "treatNorms:ei3.fNeither Agree nor Disagree" = "Norms : Neither Agree nor Disagree",
        "treatNorms:ei3.fAgree" = "Norms : Agree",
        "treatNorms:ei3.fStrongly Agree" = "Norms : Stongly Agree",
        "treatEndorsement:ei3.fStrongly Disagree" = "Endorsement : Strongly Disagree",
        "treatEndorsement:ei3.fDisagree" = "Endorsement : Disagree",
        "treatEndorsement:ei3.fNeither Agree nor Disagree" = "Endorsement : Neither Agree nor Disagree",
        "treatEndorsement:ei3.fAgree" = "Endorsement : Agree",
        "treatEndorsement:ei3.fStrongly Agree" = "Endorsement : Stongly Agree",
        "treatControl:ei3.fStrongly Disagree" = "Control : Strongly Disagree",
        "treatControl:ei3.fDisagree" = "Control : Disagree",        
        "treatControl:ei3.fNeither Agree nor Disagree" = "Control : Neither Agree nor Disagree",
        "treatControl:ei3.fAgree" = "Control : Agree",
        "treatControl:ei3.fStrongly Agree" = "Control : Stongly Agree")

modelsummary(i3,
             output = "latex",
             title = "Interaction Model: Treatment x Extrinsic / Intrinsic Motivations — Enjoys Social Aspects of Church",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Predictive Margins for EI Interaction Models
p.i1.therm <- ggpredict(i1.therm, terms = c("treat", "ei1.f"))
p.i1.p1 <- ggpredict(i1.p1, terms = c("treat", "ei1.f"))
p.i1.p2 <- ggpredict(i1.p2, terms = c("treat", "ei1.f"))
p.i1.p3 <- ggpredict(i1.p3, terms = c("treat", "ei1.f"))
p.i1.p4 <- ggpredict(i1.p4, terms = c("treat", "ei1.f"))
p.i1.behave <- ggpredict(i1.behave, terms = c("treat", "ei1.f"))
p.i1.click <- ggpredict(i1.click, terms = c("treat", "ei1.f"))

p.i1.therm$model <- "Thermometer"
p.i1.p1$model <- "Policy: Hand Washing"
p.i1.p2$model <- "Policy: Distancing"
p.i1.p3$model <- "Policy: Mask Wearing"
p.i1.p4$model <- "Policy: Prayer"
p.i1.behave$model <- "Would Click"
p.i1.click$model <- "Clicked"

p.i1.all <- rbind(p.i1.therm, p.i1.p1, p.i1.p2, p.i1.p3, p.i1.p4, p.i1.behave, p.i1.click)

p.i1.all$x <- factor(p.i1.all$x, 
                     levels = c("Norms", 
                                "Endorsement", 
                                "Control") ,
                     ordered = TRUE)

p.i1.all$model <- factor(p.i1.all$model, 
                         levels = c("Thermometer", 
                                    "Policy: Hand Washing", 
                                    "Policy: Distancing", 
                                    "Policy: Mask Wearing", 
                                    "Policy: Prayer", 
                                    "Would Click", 
                                    "Clicked"), 
                         ordered = TRUE)

p.i2.therm <- ggpredict(i2.therm, terms = c("treat", "ei2.f"))
p.i2.p1 <- ggpredict(i2.p1, terms = c("treat", "ei2.f"))
p.i2.p2 <- ggpredict(i2.p2, terms = c("treat", "ei2.f"))
p.i2.p3 <- ggpredict(i2.p3, terms = c("treat", "ei2.f"))
p.i2.p4 <- ggpredict(i2.p4, terms = c("treat", "ei2.f"))
p.i2.behave <- ggpredict(i2.behave, terms = c("treat", "ei2.f"))
p.i2.click <- ggpredict(i2.click, terms = c("treat", "ei2.f"))

p.i2.therm$model <- "Thermometer"
p.i2.p1$model <- "Policy: Hand Washing"
p.i2.p2$model <- "Policy: Distancing"
p.i2.p3$model <- "Policy: Mask Wearing"
p.i2.p4$model <- "Policy: Prayer"
p.i2.behave$model <- "Would Click"
p.i2.click$model <- "Clicked"

p.i2.all <- rbind(p.i2.therm, p.i2.p1, p.i2.p2, p.i2.p3, p.i2.p4, p.i2.behave, p.i2.click)

p.i2.all$x <- factor(p.i2.all$x, 
                     levels = c("Norms", 
                                "Endorsement", 
                                "Control") ,
                     ordered = TRUE)

p.i2.all$model <- factor(p.i2.all$model, 
                         levels = c("Thermometer", 
                                    "Policy: Hand Washing", 
                                    "Policy: Distancing", 
                                    "Policy: Mask Wearing", 
                                    "Policy: Prayer", 
                                    "Would Click", 
                                    "Clicked"), 
                         ordered = TRUE)

p.i3.therm <- ggpredict(i3.therm, terms = c("treat", "ei3.f"))
p.i3.p1 <- ggpredict(i3.p1, terms = c("treat", "ei3.f"))
p.i3.p2 <- ggpredict(i3.p2, terms = c("treat", "ei3.f"))
p.i3.p3 <- ggpredict(i3.p3, terms = c("treat", "ei3.f"))
p.i3.p4 <- ggpredict(i3.p4, terms = c("treat", "ei3.f"))
p.i3.behave <- ggpredict(i3.behave, terms = c("treat", "ei3.f"))
p.i3.click <- ggpredict(i3.click, terms = c("treat", "ei3.f"))

p.i3.therm$model <- "Thermometer"
p.i3.p1$model <- "Policy: Hand Washing"
p.i3.p2$model <- "Policy: Distancing"
p.i3.p3$model <- "Policy: Mask Wearing"
p.i3.p4$model <- "Policy: Prayer"
p.i3.behave$model <- "Would Click"
p.i3.click$model <- "Clicked"

p.i3.all <- rbind(p.i3.therm, p.i3.p1, p.i3.p2, p.i3.p3, p.i3.p4, p.i3.behave, p.i3.click)

p.i3.all$x <- factor(p.i3.all$x, 
                     levels = c("Norms", 
                                "Endorsement", 
                                "Control") ,
                     ordered = TRUE)

p.i3.all$model <- factor(p.i3.all$model, 
                         levels = c("Thermometer", 
                                    "Policy: Hand Washing", 
                                    "Policy: Distancing", 
                                    "Policy: Mask Wearing", 
                                    "Policy: Prayer", 
                                    "Would Click", 
                                    "Clicked"), 
                         ordered = TRUE)

## Margins Plots
ggplot(p.i1.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(strip.text.y = element_text(size = 7)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Extrinsic / Intrinsic Motivations - Religion Offers Comfort",
       x  = "Condition", y = "Predicted Value")

ggplot(p.i2.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(strip.text.y = element_text(size = 7)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Extrinsic / Intrinsic Motivations - Takes Religious Approach to Life",
       x  = "Condition", y = "Predicted Value")

ggplot(p.i3.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(strip.text.y = element_text(size = 7)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Extrinsic / Intrinsic Motivations - Enjoys Social Aspects of Church",
       x  = "Condition", y = "Predicted Value")

## Cleaning Up
rm(list=setdiff(ls(), c("working", "working2")))

## Religiosity Index ###########################################################
iv1 <- " ~ treat:rel1x + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"
iv2 <- " ~ treat:rel2x + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"
iv3 <- " ~ treat:rel3x + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"

## Interaction Models & Tables - Rel 1 - Regularly Attends Church
i1.therm <- lm(as.formula(paste("therm", iv1)), data = working)
i1.p1 <- lm(as.formula(paste("p1.wash", iv1)), data = working)
i1.p2 <- lm(as.formula(paste("p2.dist", iv1)), data = working)
i1.p3 <- lm(as.formula(paste("p3.mask", iv1)), data = working)
i1.p4 <- lm(as.formula(paste("p4.pray", iv1)), data = working)
i1.behave <- lm(as.formula(paste("behave", iv1)), data = working)
i1.click <- lm(as.formula(paste("clicked", iv1)), data = working2)

i1 <- list("Thermometer" = i1.therm, 
           "Policy: Hand Washing" = i1.p1, 
           "Policy: Distancing" = i1.p2, 
           "Policy: Mask Wearing" = i1.p3, 
           "Policy: Prayer" = i1.p4, 
           "Would Click" = i1.behave, 
           "Clicked" = i1.click)

cm <- c("(Intercept)" = "Intercept",
        "treatNorms:rel1xStrongly Disagree" = "Norms : Strongly Disagree",
        "treatNorms:rel1xDisagree" = "Norms : Disagree",
        "treatNorms:rel1xSomewhat Disagree" = "Norms : Somewhat Disagree",
        "treatNorms:rel1xSomewhat Agree" = "Norms : Somewhat Agree",
        "treatNorms:rel1xAgree" = "Norms : Agree",
        "treatNorms:rel1xStrongly Agree" = "Norms : Strongly Agree",
        "treatEndorsement:rel1xStrongly Disagree" = "Endorsement : Strongly Disagree",
        "treatEndorsement:rel1xDisagree" = "Endorsement : Disagree",
        "treatEndorsement:rel1xSomewhat Disagree" = "Endorsement : Somewhat Disagree",
        "treatEndorsement:rel1xSomewhat Agree" = "Endorsement : Somewhat Agree",
        "treatEndorsement:rel1xAgree" = "Endorsement : Agree",
        "treatEndorsement:rel1xStrongly Agree" = "Endorsement : Strongly Agree",
        "treatControl:rel1xStrongly Disagree" = "Control : Strongly Disagree",
        "treatControl:rel1xDisagree" = "Control : Disagree",
        "treatControl:rel1xSomewhat Disagree" = "Control : Somewhat Disagree",
        "treatControl:rel1xSomewhat Agree" = "Control : Somewhat Agree",
        "treatControl:rel1xAgree" = "Control : Agree",
        "treatControl:rel1xStrongly Agree" = "Control : Strongly Agree")

modelsummary(i1,
             output = "latex",
             title = "Treatment x Regularly Attends Church",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Interaction Models & Tables - Rel 2 - Spiritual Values
i2.therm <- lm(as.formula(paste("therm", iv2)), data = working)
i2.p1 <- lm(as.formula(paste("p1.wash", iv2)), data = working)
i2.p2 <- lm(as.formula(paste("p2.dist", iv2)), data = working)
i2.p3 <- lm(as.formula(paste("p3.mask", iv2)), data = working)
i2.p4 <- lm(as.formula(paste("p4.pray", iv2)), data = working)
i2.behave <- lm(as.formula(paste("behave", iv2)), data = working)
i2.click <- lm(as.formula(paste("clicked", iv2)), data = working2)

i2 <- list("Thermometer" = i2.therm, 
           "Policy: Hand Washing" = i2.p1, 
           "Policy: Distancing" = i2.p2, 
           "Policy: Mask Wearing" = i2.p3, 
           "Policy: Prayer" = i2.p4, 
           "Would Click" = i2.behave, 
           "Clicked" = i2.click)

cm <- c("(Intercept)" = "Intercept",
        "treatNorms:rel2xStrongly Disagree" = "Norms : Strongly Disagree",
        "treatNorms:rel2xDisagree" = "Norms : Disagree",
        "treatNorms:rel2xSomewhat Disagree" = "Norms : Somewhat Disagree",
        "treatNorms:rel2xSomewhat Agree" = "Norms : Somewhat Agree",
        "treatNorms:rel2xAgree" = "Norms : Agree",
        "treatNorms:rel2xStrongly Agree" = "Norms : Strongly Agree",
        "treatEndorsement:rel2xStrongly Disagree" = "Endorsement : Strongly Disagree",
        "treatEndorsement:rel2xDisagree" = "Endorsement : Disagree",
        "treatEndorsement:rel2xSomewhat Disagree" = "Endorsement : Somewhat Disagree",
        "treatEndorsement:rel2xSomewhat Agree" = "Endorsement : Somewhat Agree",
        "treatEndorsement:rel2xAgree" = "Endorsement : Agree",
        "treatEndorsement:rel2xStrongly Agree" = "Endorsement : Strongly Agree",
        "treatControl:rel2xStrongly Disagree" = "Control : Strongly Disagree",
        "treatControl:rel2xDisagree" = "Control : Disagree",
        "treatControl:rel2xSomewhat Disagree" = "Control : Somewhat Disagree",
        "treatControl:rel2xSomewhat Agree" = "Control : Somewhat Agree",
        "treatControl:rel2xAgree" = "Control : Agree",
        "treatControl:rel2xStrongly Agree" = "Control : Strongly Agree")

modelsummary(i2,
             output = "latex",
             title = "Treatment x Spiritual Values",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Interaction Models & Tables - Rel 3 - Americans Should be More Religious
i3.therm <- lm(as.formula(paste("therm", iv3)), data = working)
i3.p1 <- lm(as.formula(paste("p1.wash", iv3)), data = working)
i3.p2 <- lm(as.formula(paste("p2.dist", iv3)), data = working)
i3.p3 <- lm(as.formula(paste("p3.mask", iv3)), data = working)
i3.p4 <- lm(as.formula(paste("p4.pray", iv3)), data = working)
i3.behave <- lm(as.formula(paste("behave", iv3)), data = working)
i3.click <- lm(as.formula(paste("clicked", iv3)), data = working2)

i3 <- list("Thermometer" = i3.therm, 
           "Policy: Hand Washing" = i3.p1, 
           "Policy: Distancing" = i3.p2, 
           "Policy: Mask Wearing" = i3.p3, 
           "Policy: Prayer" = i3.p4, 
           "Would Click" = i3.behave, 
           "Clicked" = i3.click)

cm <- c("(Intercept)" = "Intercept",
        "treatNorms:rel3xStrongly Disagree" = "Norms : Strongly Disagree",
        "treatNorms:rel3xDisagree" = "Norms : Disagree",
        "treatNorms:rel3xSomewhat Disagree" = "Norms : Somewhat Disagree",
        "treatNorms:rel3xSomewhat Agree" = "Norms : Somewhat Agree",
        "treatNorms:rel3xAgree" = "Norms : Agree",
        "treatNorms:rel3xStrongly Agree" = "Norms : Strongly Agree",
        "treatEndorsement:rel3xStrongly Disagree" = "Endorsement : Strongly Disagree",
        "treatEndorsement:rel3xDisagree" = "Endorsement : Disagree",
        "treatEndorsement:rel3xSomewhat Disagree" = "Endorsement : Somewhat Disagree",
        "treatEndorsement:rel3xSomewhat Agree" = "Endorsement : Somewhat Agree",
        "treatEndorsement:rel3xAgree" = "Endorsement : Agree",
        "treatEndorsement:rel3xStrongly Agree" = "Endorsement : Strongly Agree",
        "treatControl:rel3xStrongly Disagree" = "Control : Strongly Disagree",
        "treatControl:rel3xDisagree" = "Control : Disagree",
        "treatControl:rel3xSomewhat Disagree" = "Control : Somewhat Disagree",
        "treatControl:rel3xSomewhat Agree" = "Control : Somewhat Agree",
        "treatControl:rel3xAgree" = "Control : Agree",
        "treatControl:rel3xStrongly Agree" = "Control : Strongly Agree") 

modelsummary(i3,
             output = "latex",
             title = "Treatment × Americans Should be More Religious",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Predictive Margins for Religiosity Index Interaction Models
p.i1.therm <- ggpredict(i1.therm, terms = c("treat", "rel1x"))
p.i1.p1 <- ggpredict(i1.p1, terms = c("treat", "rel1x"))
p.i1.p2 <- ggpredict(i1.p2, terms = c("treat", "rel1x"))
p.i1.p3 <- ggpredict(i1.p3, terms = c("treat", "rel1x"))
p.i1.p4 <- ggpredict(i1.p4, terms = c("treat", "rel1x"))
p.i1.behave <- ggpredict(i1.behave, terms = c("treat", "rel1x"))
p.i1.click <- ggpredict(i1.click, terms = c("treat", "rel1x"))

p.i1.therm$model <- "Thermometer"
p.i1.p1$model <- "Policy: Hand Washing"
p.i1.p2$model <- "Policy: Distancing"
p.i1.p3$model <- "Policy: Mask Wearing"
p.i1.p4$model <- "Policy: Prayer"
p.i1.behave$model <- "Would Click"
p.i1.click$model <- "Clicked"

p.i1.all <- rbind(p.i1.therm, p.i1.p1, p.i1.p2, p.i1.p3, p.i1.p4, p.i1.behave, p.i1.click)

p.i1.all$x <- factor(p.i1.all$x, 
                     levels = c("Norms", 
                                "Endorsement", 
                                "Control") ,
                     ordered = TRUE)

p.i1.all$model <- factor(p.i1.all$model, 
                         levels = c("Thermometer", 
                                    "Policy: Hand Washing", 
                                    "Policy: Distancing", 
                                    "Policy: Mask Wearing", 
                                    "Policy: Prayer", 
                                    "Would Click", 
                                    "Clicked"), 
                         ordered = TRUE)

p.i2.therm <- ggpredict(i2.therm, terms = c("treat", "rel2x"))
p.i2.p1 <- ggpredict(i2.p1, terms = c("treat", "rel2x"))
p.i2.p2 <- ggpredict(i2.p2, terms = c("treat", "rel2x"))
p.i2.p3 <- ggpredict(i2.p3, terms = c("treat", "rel2x"))
p.i2.p4 <- ggpredict(i2.p4, terms = c("treat", "rel2x"))
p.i2.behave <- ggpredict(i2.behave, terms = c("treat", "rel2x"))
p.i2.click <- ggpredict(i2.click, terms = c("treat", "rel2x"))

p.i2.therm$model <- "Thermometer"
p.i2.p1$model <- "Policy: Hand Washing"
p.i2.p2$model <- "Policy: Distancing"
p.i2.p3$model <- "Policy: Mask Wearing"
p.i2.p4$model <- "Policy: Prayer"
p.i2.behave$model <- "Would Click"
p.i2.click$model <- "Clicked"

p.i2.all <- rbind(p.i2.therm, p.i2.p1, p.i2.p2, p.i2.p3, p.i2.p4, p.i2.behave, p.i2.click)

p.i2.all$x <- factor(p.i2.all$x, 
                     levels = c("Norms", 
                                "Endorsement", 
                                "Control") ,
                     ordered = TRUE)

p.i2.all$model <- factor(p.i2.all$model, 
                         levels = c("Thermometer", 
                                    "Policy: Hand Washing", 
                                    "Policy: Distancing", 
                                    "Policy: Mask Wearing", 
                                    "Policy: Prayer", 
                                    "Would Click", 
                                    "Clicked"), 
                         ordered = TRUE)

p.i3.therm <- ggpredict(i3.therm, terms = c("treat", "rel3x"))
p.i3.p1 <- ggpredict(i3.p1, terms = c("treat", "rel3x"))
p.i3.p2 <- ggpredict(i3.p2, terms = c("treat", "rel3x"))
p.i3.p3 <- ggpredict(i3.p3, terms = c("treat", "rel3x"))
p.i3.p4 <- ggpredict(i3.p4, terms = c("treat", "rel3x"))
p.i3.behave <- ggpredict(i3.behave, terms = c("treat", "rel3x"))
p.i3.click <- ggpredict(i3.click, terms = c("treat", "rel3x"))

p.i3.therm$model <- "Thermometer"
p.i3.p1$model <- "Policy: Hand Washing"
p.i3.p2$model <- "Policy: Distancing"
p.i3.p3$model <- "Policy: Mask Wearing"
p.i3.p4$model <- "Policy: Prayer"
p.i3.behave$model <- "Would Click"
p.i3.click$model <- "Clicked"

p.i3.all <- rbind(p.i3.therm, p.i3.p1, p.i3.p2, p.i3.p3, p.i3.p4, p.i3.behave, p.i3.click)

p.i3.all$x <- factor(p.i3.all$x, 
                     levels = c("Norms", 
                                "Endorsement", 
                                "Control") ,
                     ordered = TRUE)

p.i3.all$model <- factor(p.i3.all$model, 
                         levels = c("Thermometer", 
                                    "Policy: Hand Washing", 
                                    "Policy: Distancing", 
                                    "Policy: Mask Wearing", 
                                    "Policy: Prayer", 
                                    "Would Click", 
                                    "Clicked"), 
                         ordered = TRUE)


## Plots
ggplot(p.i1.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(strip.text.y = element_text(size = 7)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Regularly Attends Church",
       x  = "Condition", y = "Predicted Value")

ggplot(p.i2.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(strip.text.y = element_text(size = 7)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Spiritual Values",
       x  = "Condition", y = "Predicted Value")

ggplot(p.i3.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(strip.text.y = element_text(size = 7)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Americans Should be More Religious",
       x  = "Condition", y = "Predicted Value")

## Cleaning Up
rm(list=setdiff(ls(), "working"))


# Exploratory Analysis of Partisanship #########################################
## Conditional Average Treatment Effects #######################################
levels(working$party.f)   # Change the reference category to independent
working$party.f2 <- factor(working$party.f, levels = c("Strong Democrat", "Democrat", 
                                                       "Lean Democrat", "Independent",
                                                       "Lean Republican", 
                                                       "Republican", "Strong Republican"))
working$party.32 <- factor(working$party.3, levels = c("Democrat",
                                                       "Independent", 
                                                       "Republican"))

w2 <- working %>% 
  filter(behave > 4)

iv1 <- " ~ treat * party.32 + gender + age + education + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"
iv2 <- " ~ treat * party.f2 + gender + age + education + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"


### Three Category Party Analysis ##############################################
## Simplified Model, No Controls
# a1 <- lm(therm ~ treat * party.32, data = working)
# a2 <- lm(p1.wash ~ treat * party.32, data = working)
# a3 <- lm(p2.dist ~ treat * party.32, data = working)
# a4 <- lm(p3.mask ~ treat * party.32, data = working)
# a5 <- lm(p4.pray ~ treat * party.32, data = working)
# a6 <- lm(behave ~ treat * party.32, data = working)
# a7 <- lm(clicked ~ treat * party.32, data = w2)

## Model With All Controls 
a1 <- lm(as.formula(paste("therm", iv1)), data = working)
a2 <- lm(as.formula(paste("p1.wash", iv1)), data = working)
a3 <- lm(as.formula(paste("p2.dist", iv1)), data = working)
a4 <- lm(as.formula(paste("p3.mask", iv1)), data = working)
a5 <- lm(as.formula(paste("p4.pray", iv1)), data = working)
a6 <- lm(as.formula(paste("behave", iv1)), data = working)
a7 <- lm(as.formula(paste("clicked", iv1)), data = w2)

### 3 Cat Margins & Margins Plots
int.therm <- ggpredict(a1, terms = c("treat", "party.32"))
int.p1 <- ggpredict(a2, terms = c("treat", "party.32"))
int.p2 <- ggpredict(a3, terms = c("treat", "party.32"))
int.p3 <- ggpredict(a4, terms = c("treat", "party.32"))
int.p4 <- ggpredict(a5, terms = c("treat", "party.32"))
int.behave <- ggpredict(a6, terms = c("treat", "party.32"))
int.click <- ggpredict(a7, terms = c("treat", "party.32"))

int.therm$model <- "Thermometer"
int.p1$model <- "Policy: Hand Washing"
int.p2$model <- "Policy: Distancing"
int.p3$model <- "Policy: Mask Wearing"
int.p4$model <- "Policy: Prayer"
int.behave$model <- "Would Click"
int.click$model <- "Clicked"

int.all <- rbind(int.therm, int.p1, int.p2, int.p3, int.p4, int.behave, int.click)

int.all$model <- factor(int.all$model, 
                        levels = c("Thermometer", 
                                   "Policy: Hand Washing", 
                                   "Policy: Distancing", 
                                   "Policy: Mask Wearing", 
                                   "Policy: Prayer", 
                                   "Would Click", 
                                   "Clicked"), 
                        ordered = TRUE)

ggplot(int.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model),
             scales = "free_y") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Party (Three Category)",
       x  = "Condition", y = "Predicted Value",
       caption = "Interaction model with controls.")

rm(int.all, int.therm, int.p1, int.p2, int.p3, int.p4, int.behave, int.click)

## 3 Cat Contrasts & PWPPs 
## Use this for models without controls:
# z.a1 <- emmeans(a1, ~ treat*party.32)
# z.a2 <- emmeans(a2, ~ treat*party.32)
# z.a3 <- emmeans(a3, ~ treat*party.32)
# z.a4 <- emmeans(a4, ~ treat*party.32)
# z.a5 <- emmeans(a5, ~ treat*party.32)
# z.a6 <- emmeans(a6, ~ treat*party.32)
# z.a7 <- emmeans(a7, ~ treat*party.32)

## Use this for models with controls:
z.a1 <- emmeans(a1, ~ treat*party.32, non.nuis = c("treat", "party.32"))
z.a2 <- emmeans(a2, ~ treat*party.32, non.nuis = c("treat", "party.32"))
z.a3 <- emmeans(a3, ~ treat*party.32, non.nuis = c("treat", "party.32"))
z.a4 <- emmeans(a4, ~ treat*party.32, non.nuis = c("treat", "party.32"))
z.a5 <- emmeans(a5, ~ treat*party.32, non.nuis = c("treat", "party.32"))
z.a6 <- emmeans(a6, ~ treat*party.32, non.nuis = c("treat", "party.32"))
z.a7 <- emmeans(a7, ~ treat*party.32, non.nuis = c("treat", "party.32"))

pwpp(z.a1, by = "party.32", 
     xlab = "Tukey-Adjusted P Value",
     ylab = "Estimated Marginal Mean") +
  ggplot2::theme_bw() +
  ggplot2::ggtitle("Pairwise Comparisons", 
                   subtitle = "Thermometer") +
  ggplot2::facet_grid(~ party.32)

pwpp(z.a2, by = "party.32", 
     xlab = "Tukey-Adjusted P Value",
     ylab = "Estimated Marginal Mean") +
  ggplot2::theme_bw() +
  ggplot2::ggtitle("Pairwise Comparisons", 
                   subtitle = "Policy: Hand Washing") +
  ggplot2::facet_grid(~ party.32)

pwpp(z.a3, by = "party.32", 
     xlab = "Tukey-Adjusted P Value",
     ylab = "Estimated Marginal Mean") +
  ggplot2::theme_bw() +
  ggplot2::ggtitle("Pairwise Comparisons", 
                   subtitle = "Policy: Distancing") +
  ggplot2::facet_grid(~ party.32)

pwpp(z.a4, by = "party.32", 
     xlab = "Tukey-Adjusted P Value",
     ylab = "Estimated Marginal Mean") +
  ggplot2::theme_bw() +
  ggplot2::ggtitle("Pairwise Comparisons", 
                   subtitle = "Policy: Mask Wearing") +
  ggplot2::facet_grid(~ party.32)

pwpp(z.a5, by = "party.32", 
     xlab = "Tukey-Adjusted P Value",
     ylab = "Estimated Marginal Mean") +
  ggplot2::theme_bw() +
  ggplot2::ggtitle("Pairwise Comparisons", 
                   subtitle = "Policy: Prayer") +
  ggplot2::facet_grid(~ party.32)

pwpp(z.a6, by = "party.32", 
     xlab = "Tukey-Adjusted P Value",
     ylab = "Estimated Marginal Mean") +
  ggplot2::theme_bw() +
  ggplot2::ggtitle("Pairwise Comparisons", 
                   subtitle = "Would Click") +
  ggplot2::facet_grid(~ party.32)

pwpp(z.a7, by = "party.32", 
     xlab = "Tukey-Adjusted P Value",
     ylab = "Estimated Marginal Mean") +
  ggplot2::theme_bw() +
  ggplot2::ggtitle("Pairwise Comparisons", 
                   subtitle = "Clicked") +
  ggplot2::facet_grid(~ party.32)


### Seven Category Party Analysis ##############################################
## Simplified Model, No Controls
# a1 <- lm(therm ~ treat * party.f2, data = working)
# a2 <- lm(p1.wash ~ treat * party.f2, data = working)
# a3 <- lm(p2.dist ~ treat * party.f2, data = working)
# a4 <- lm(p3.mask ~ treat * party.f2, data = working)
# a5 <- lm(p4.pray ~ treat * party.f2, data = working)
# a6 <- lm(behave ~ treat * party.f2, data = working)
# a7 <- lm(clicked ~ treat * party.f2, data = w2)

## Model With All Controls
a1 <- lm(as.formula(paste("therm", iv2)), data = working)
a2 <- lm(as.formula(paste("p1.wash", iv2)), data = working)
a3 <- lm(as.formula(paste("p2.dist", iv2)), data = working)
a4 <- lm(as.formula(paste("p3.mask", iv2)), data = working)
a5 <- lm(as.formula(paste("p4.pray", iv2)), data = working)
a6 <- lm(as.formula(paste("behave", iv2)), data = working)
a7 <- lm(as.formula(paste("clicked", iv2)), data = w2)

## 7 Cat Margins & Margins Plots
int.therm <- ggpredict(a1, terms = c("treat", "party.f2"))
int.p1 <- ggpredict(a2, terms = c("treat", "party.f2"))
int.p2 <- ggpredict(a3, terms = c("treat", "party.f2"))
int.p3 <- ggpredict(a4, terms = c("treat", "party.f2"))
int.p4 <- ggpredict(a5, terms = c("treat", "party.f2"))
int.behave <- ggpredict(a6, terms = c("treat", "party.f2"))
int.click <- ggpredict(a7, terms = c("treat", "party.f2"))

int.therm$model <- "Thermometer"
int.p1$model <- "Policy: Hand Washing"
int.p2$model <- "Policy: Distancing"
int.p3$model <- "Policy: Mask Wearing"
int.p4$model <- "Policy: Prayer"
int.behave$model <- "Would Click"
int.click$model <- "Clicked"

int.all <- rbind(int.therm, int.p1, int.p2, int.p3, int.p4, int.behave, int.click)

int.all$model <- factor(int.all$model, 
                        levels = c("Thermometer", 
                                   "Policy: Hand Washing", 
                                   "Policy: Distancing", 
                                   "Policy: Mask Wearing", 
                                   "Policy: Prayer", 
                                   "Would Click", 
                                   "Clicked"), 
                        ordered = TRUE)

ggplot(int.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model),
             scales = "free_y") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Party (Seven Category)",
       x  = "Condition", y = "Predicted Value",
       caption = "Interaction model with controls.") +
  theme(strip.text.y = element_text(size = 7))

## Cleaning Up
rm(list=setdiff(ls(), "working"))

## Subgroup Analysis ###########################################################
w <- working %>% 
  filter(party.3 == "Republican")

w2 <- working %>% 
  filter(behave > 4)

iv1 <- " ~ treat + gender + age + education + party.f + covid.threat + covid.tested + covid.friends"
iv2 <- " ~ treat + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends"


### Baseline OLS Models ########################################################
## The Models
a1 <- lm(therm ~ treat, data = w)
a2 <- lm(p1.wash ~ treat, data = w)
a3 <- lm(p2.dist ~ treat, data = w)
a4 <- lm(p3.mask ~ treat, data = w)
a5 <- lm(p4.pray ~ treat, data = w)
a6 <- lm(behave ~ treat, data = working)
a7 <- lm(clicked ~ treat, data = w2)

## Create Dataframe for Coefficient Plot
x1 <- tidy(a1, conf.int = TRUE)
x2 <- tidy(a2, conf.int = TRUE)
x3 <- tidy(a3, conf.int = TRUE)
x4 <- tidy(a4, conf.int = TRUE)
x5 <- tidy(a5, conf.int = TRUE)
x6 <- tidy(a6, conf.int = TRUE)
x7 <- tidy(a7, conf.int = TRUE)

x1$out <- "Thermometer"
x2$out <- "Policy: Hand Washing"
x3$out <- "Policy: Distancing"
x4$out <- "Policy: Mask Wearing"
x5$out <- "Policy: Prayer"
x6$out <- "Would Click"
x7$out <- "Clicked"

z1 <- tidy(a1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(a2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(a3, conf.int = TRUE, conf.level = .9)
z4 <- tidy(a4, conf.int = TRUE, conf.level = .9)
z5 <- tidy(a5, conf.int = TRUE, conf.level = .9)
z6 <- tidy(a6, conf.int = TRUE, conf.level = .9)
z7 <- tidy(a7, conf.int = TRUE, conf.level = .9)

z1$out <- "Thermometer"
z2$out <- "Policy: Hand Washing"
z3$out <- "Policy: Distancing"
z4$out <- "Policy: Mask Wearing"
z5$out <- "Policy: Prayer"
z6$out <- "Would Click"
z7$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)
z4 <- z4 %>% rename(low = conf.low, high = conf.high)
z5 <- z5 %>% rename(low = conf.low, high = conf.high)
z6 <- z6 %>% rename(low = conf.low, high = conf.high)
z7 <- z7 %>% rename(low = conf.low, high = conf.high)

## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3, z4, z5, z6, z7)
x <- rbind(x1, x2, x3, x4, x5, x6, x7)

q <- z %>% full_join(x)

q$out <- factor(q$out, 
                levels = c("Thermometer", 
                           "Policy: Hand Washing", 
                           "Policy: Distancing", 
                           "Policy: Mask Wearing", 
                           "Policy: Prayer", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

## Filter out Intercept
y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
ggplot(y) +
  geom_hline(yintercept=0, lty="11", colour="grey30") +
  geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.2) +
  geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
  geom_point(aes(cond, estimate)) +
  labs(colour="", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out)) + 
  labs(title = "Coefficient Plot: Baseline OLS Model",
       caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) 

y1 <- y
y1$mdl <- "Baseline"


### Demographic OLS Models #####################################################
rm(q, x, x1, x2, x3, x4, x5, x6, x7, y, z, z1, z2, z3, z4, z5, z6, z7)

## The Models
b1 <- lm(as.formula(paste("therm", iv1)), data = w)
b2 <- lm(as.formula(paste("p1.wash", iv1)), data = w)
b3 <- lm(as.formula(paste("p2.dist", iv1)), data = w)
b4 <- lm(as.formula(paste("p3.mask", iv1)), data = w)
b5 <- lm(as.formula(paste("p4.pray", iv1)), data = w)
b6 <- lm(as.formula(paste("behave", iv1)), data = w)
b7 <- lm(as.formula(paste("clicked", iv1)), data = w2)

## Create Dataframe for Coefficient Plot
x1 <- tidy(b1, conf.int = TRUE)
x2 <- tidy(b2, conf.int = TRUE)
x3 <- tidy(b3, conf.int = TRUE)
x4 <- tidy(b4, conf.int = TRUE)
x5 <- tidy(b5, conf.int = TRUE)
x6 <- tidy(b6, conf.int = TRUE)
x7 <- tidy(b7, conf.int = TRUE)

x1$out <- "Thermometer"
x2$out <- "Policy: Hand Washing"
x3$out <- "Policy: Distancing"
x4$out <- "Policy: Mask Wearing"
x5$out <- "Policy: Prayer"
x6$out <- "Would Click"
x7$out <- "Clicked"

z1 <- tidy(b1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(b2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(b3, conf.int = TRUE, conf.level = .9)
z4 <- tidy(b4, conf.int = TRUE, conf.level = .9)
z5 <- tidy(b5, conf.int = TRUE, conf.level = .9)
z6 <- tidy(b6, conf.int = TRUE, conf.level = .9)
z7 <- tidy(b7, conf.int = TRUE, conf.level = .9)

z1$out <- "Thermometer"
z2$out <- "Policy: Hand Washing"
z3$out <- "Policy: Distancing"
z4$out <- "Policy: Mask Wearing"
z5$out <- "Policy: Prayer"
z6$out <- "Would Click"
z7$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)
z4 <- z4 %>% rename(low = conf.low, high = conf.high)
z5 <- z5 %>% rename(low = conf.low, high = conf.high)
z6 <- z6 %>% rename(low = conf.low, high = conf.high)
z7 <- z7 %>% rename(low = conf.low, high = conf.high)

## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3, z4, z5, z6, z7)
x <- rbind(x1, x2, x3, x4, x5, x6, x7)

q <- z %>% full_join(x)

q <- q %>%
  filter(term == "(Intercept)" | term == "treatEndorsement" | term == "treatNorms")

q$out <- factor(q$out, 
                levels = c("Thermometer", 
                           "Policy: Hand Washing", 
                           "Policy: Distancing", 
                           "Policy: Mask Wearing", 
                           "Policy: Prayer", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

## Filter out Intercept
y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
ggplot(y) +
  geom_hline(yintercept=0, lty="11", colour="grey30") +
  geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.2) +
  geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
  geom_point(aes(cond, estimate)) +
  labs(colour="", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out)) + 
  labs(title = "Coefficient Plot: OLS Model with Demographic Controls",
       caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) 

y2 <- y
y2$mdl <- "Demographic Controls"


### OLS Models All Controls ####################################################
rm(q, x, x1, x2, x3, x4, x5, x6, x7, y, z, z1, z2, z3, z4, z5, z6, z7)

## The Models
c1 <- lm(as.formula(paste("therm", iv2)), data = w)
c2 <- lm(as.formula(paste("p1.wash", iv2)), data = w)
c3 <- lm(as.formula(paste("p2.dist", iv2)), data = w)
c4 <- lm(as.formula(paste("p3.mask", iv2)), data = w)
c5 <- lm(as.formula(paste("p4.pray", iv2)), data = w)
c6 <- lm(as.formula(paste("behave", iv2)), data = w)
c7 <- lm(as.formula(paste("clicked", iv2)), data = w2)

## Create Dataframe for Coefficient Plot
x1 <- tidy(c1, conf.int = TRUE)
x2 <- tidy(c2, conf.int = TRUE)
x3 <- tidy(c3, conf.int = TRUE)
x4 <- tidy(c4, conf.int = TRUE)
x5 <- tidy(c5, conf.int = TRUE)
x6 <- tidy(c6, conf.int = TRUE)
x7 <- tidy(c7, conf.int = TRUE)

x1$out <- "Thermometer"
x2$out <- "Policy: Hand Washing"
x3$out <- "Policy: Distancing"
x4$out <- "Policy: Mask Wearing"
x5$out <- "Policy: Prayer"
x6$out <- "Would Click"
x7$out <- "Clicked"

z1 <- tidy(c1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(c2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(c3, conf.int = TRUE, conf.level = .9)
z4 <- tidy(c4, conf.int = TRUE, conf.level = .9)
z5 <- tidy(c5, conf.int = TRUE, conf.level = .9)
z6 <- tidy(c6, conf.int = TRUE, conf.level = .9)
z7 <- tidy(c7, conf.int = TRUE, conf.level = .9)

z1$out <- "Thermometer"
z2$out <- "Policy: Hand Washing"
z3$out <- "Policy: Distancing"
z4$out <- "Policy: Mask Wearing"
z5$out <- "Policy: Prayer"
z6$out <- "Would Click"
z7$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)
z4 <- z4 %>% rename(low = conf.low, high = conf.high)
z5 <- z5 %>% rename(low = conf.low, high = conf.high)
z6 <- z6 %>% rename(low = conf.low, high = conf.high)
z7 <- z7 %>% rename(low = conf.low, high = conf.high)

## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3, z4, z5, z6, z7)
x <- rbind(x1, x2, x3, x4, x5, x6, x7)

q <- z %>% full_join(x)

q <- q %>%
  filter(term == "(Intercept)" | term == "treatEndorsement" | term == "treatNorms")

q$out <- factor(q$out, 
                levels = c("Thermometer", 
                           "Policy: Hand Washing", 
                           "Policy: Distancing", 
                           "Policy: Mask Wearing", 
                           "Policy: Prayer", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

## Filter out Intercept
y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
ggplot(y) +
  geom_hline(yintercept=0, lty="11", colour="grey30") +
  geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.2) +
  geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
  geom_point(aes(cond, estimate)) +
  labs(colour="", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out)) + 
  labs(title = "Coefficient Plot: OLS Model with Controls - Republican Sub-Sample",
       caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

y3 <- y
y3$mdl <- "All Controls"


### Grid Plot of All Models ####################################################
rm(q, x, x1, x2, x3, x4, x5, x6, x7, y, z, z1, z2, z3, z4, z5, z6, z7)

y <- rbind(y1, y2, y3)

y$mdl <- factor(y$mdl, 
                levels = c("Baseline", 
                           "Demographic Controls", 
                           "All Controls"), 
                ordered = TRUE)

ggplot(y) +
  geom_hline(yintercept=0, lty="11", colour="grey30") +
  geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.2) +
  geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
  geom_point(aes(cond, estimate)) +
  labs(colour="", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out), rows = vars(mdl)) + 
  labs(title = "Coefficient Plot: OLS Models - Republican Sub-Sample",
       caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

## Cleaning Up
rm(list=setdiff(ls(), "working"))


# Synthetic Index of Attitudinal Outcomes ######################################
## PCA and Scores ##############################################################
pc1 <- princomp(~ therm + p1.wash + p2.dist + p3.mask + p4.pray, 
                data = working, na.action = na.exclude, cor = TRUE)

summary(pc1)
working$pc.score <- pc1$scores[,1] 


## New OLS with PCA as Outcome #################################################
iv1 <- " ~ treat + gender + age + education + party.f + covid.threat + covid.tested + covid.friends"
iv2 <- " ~ treat + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends"

w2 <- working %>% 
  filter(behave > 4)

### Baseline OLS Models ########################################################
## The Models
a1 <- lm(pc.score ~ treat, data = working)
a2 <- lm(behave ~ treat, data = working)
a3 <- lm(clicked ~ treat, data = w2)

## Create Dataframe for Coefficient Plot
x1 <- tidy(a1, conf.int = TRUE)
x2 <- tidy(a2, conf.int = TRUE)
x3 <- tidy(a3, conf.int = TRUE)

x1$out <- "Attitudinal Outcomes"
x2$out <- "Would Click"
x3$out <- "Clicked"

z1 <- tidy(a1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(a2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(a3, conf.int = TRUE, conf.level = .9)

z1$out <- "Attitudinal Outcomes"
z2$out <- "Would Click"
z3$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)

## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3)
x <- rbind(x1, x2, x3)

q <- z %>% full_join(x)

q$out <- factor(q$out, 
                levels = c("Attitudinal Outcomes", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

## Filter out Intercept
y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
ggplot(y) +
  geom_hline(yintercept=0, lty="11", colour="grey30") +
  geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.2) +
  geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
  geom_point(aes(cond, estimate)) +
  labs(colour="", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out)) + 
  labs(title = "Coefficient Plot: Baseline OLS Model",
       caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  scale_y_continuous(breaks = c(-0.2, 0, 0.2, 0.4, 0.6, 0.8), limits = c(-0.2,0.8))

y1 <- y
y1$mdl <- "Baseline"

### Demographic OLS Models #####################################################
rm(q, x, x1, x2, x3, y, z, z1, z2, z3)

## The Models
b1 <- lm(as.formula(paste("pc.score", iv1)), data = working)
b2 <- lm(as.formula(paste("behave", iv1)), data = working)
b3 <- lm(as.formula(paste("clicked", iv1)), data = working)


## Create Dataframe for Coefficient Plot
x1 <- tidy(b1, conf.int = TRUE)
x2 <- tidy(b2, conf.int = TRUE)
x3 <- tidy(b3, conf.int = TRUE)

x1$out <- "Attitudinal Outcomes"
x2$out <- "Would Click"
x3$out <- "Clicked"

z1 <- tidy(b1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(b2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(b3, conf.int = TRUE, conf.level = .9)

z1$out <- "Attitudinal Outcomes"
z2$out <- "Would Click"
z3$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)

## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3)
x <- rbind(x1, x2, x3)

q <- z %>% full_join(x)

q <- q %>%
  filter(term == "(Intercept)" | term == "treatEndorsement" | term == "treatNorms")

q$out <- factor(q$out, 
                levels = c("Attitudinal Outcomes", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

## Filter out Intercept
y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
ggplot(y) +
  geom_hline(yintercept=0, lty="11", colour="grey30") +
  geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.2) +
  geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
  geom_point(aes(cond, estimate)) +
  labs(colour="", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out)) + 
  labs(title = "Coefficient Plot: OLS Model with Demographic Controls",
       caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  scale_y_continuous(breaks = c(-0.2, 0, 0.2, 0.4, 0.6, 0.8), limits = c(-0.2,0.8))

y2 <- y
y2$mdl <- "Demographic Controls"

### OLS Models All Controls ####################################################
rm(q, x, x1, x2, x3, y, z, z1, z2, z3)

## The Models
c1 <- lm(as.formula(paste("pc.score", iv2)), data = working)
c2 <- lm(as.formula(paste("behave", iv2)), data = working)
c3 <- lm(as.formula(paste("clicked", iv2)), data = working)

## Create Dataframe for Coefficient Plot
x1 <- tidy(c1, conf.int = TRUE)
x2 <- tidy(c2, conf.int = TRUE)
x3 <- tidy(c3, conf.int = TRUE)

x1$out <- "Attitudinal Outcomes"
x2$out <- "Would Click"
x3$out <- "Clicked"

z1 <- tidy(c1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(c2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(c3, conf.int = TRUE, conf.level = .9)

z1$out <- "Attitudinal Outcomes"
z2$out <- "Would Click"
z3$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)


## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3)
x <- rbind(x1, x2, x3)

q <- z %>% full_join(x)

q <- q %>%
  filter(term == "(Intercept)" | term == "treatEndorsement" | term == "treatNorms")

q$out <- factor(q$out, 
                levels = c("Attitudinal Outcomes", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

## Filter out Intercept
y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
ggplot(y) +
  geom_hline(yintercept=0, lty="11", colour="grey30") +
  geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.2) +
  geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
  geom_point(aes(cond, estimate)) +
  labs(colour="", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out)) + 
  labs(title = "Coefficient Plot: OLS Model with Controls",
       caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  scale_y_continuous(breaks = c(-0.2, 0, 0.2, 0.4, 0.6, 0.8), limits = c(-0.2,0.8))

y3 <- y
y3$mdl <- "All Controls"

### Grid Plot of All Models ####################################################
rm(q, x, x1, x2, x3, y, z, z1, z2, z3)

y <- rbind(y1, y2, y3)

y$mdl <- factor(y$mdl, 
                levels = c("Baseline", 
                           "Demographic Controls", 
                           "All Controls"), 
                ordered = TRUE)

ggplot(y) +
  geom_hline(yintercept=0, lty="11", colour="grey30") +
  geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.075) +
  geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
  geom_point(aes(cond, estimate)) +
  labs(colour="", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out), rows = vars(mdl)) + 
  labs(title = "Coefficient Plot: OLS Models",
       caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

## Cleaning Up
rm(list=setdiff(ls(), "working"))


# Assessing Anchoring Effects ##################################################
pols3 <- lm(p3.mask ~ p1.wash, data = working)
summary(pols3)

pols4 <- lm(p3.mask ~ p1.wash + p2.dist, data = working)
summary(pols3)

pols5 <- lm(p3.mask ~ p1.wash + p2.dist + treat, data = working)
summary(pols5)

xmodels <- list(pols3, pols4, pols5)

cm = c("(Intercept)" = "Intercept",
       "p1.wash" = "Policy: Hand Washing (Item 1)",
       "p2.dist" = "Policy: Distancing (Item 2)",
       "treatEndorsement" = "Endorsement",
       "treatNorms" = "Norms")

modelsummary(xmodels, 
             title = 'Mask Wearing (Item 3)',
             coef_map = cm,
             stars = TRUE, 
             gof_omit = 'IC|Log|Adj')

# Prior Knowledge ##############################################################
rm(list=setdiff(ls(), "working"))

working2 <- working %>% 
  unite("knew", 13:15, na.rm = TRUE, remove = FALSE)

working2$knew <- as.numeric(working2$knew)

working2$knewit <- dplyr::recode(working2$knew,
                                 `3` = "Knew It",
                                 `2` = "Didn't Know It",
                                 `1` = "Don't Believe It",
                                 .default = NA_character_)

working2$knewit <- as.factor(working2$knewit)

t <- table(working2$treat, working2$knewit)
t
prop.table(t, 1)

ggplot(data = working2, aes(x = factor(knewit))) +
  geom_bar(stat = "count") +
  geom_text(stat = "count", aes(label = ..count..), vjust = -0.5, size = 2.5) +
  labs(title = "Prior Knowledge of Treatment Information, by Treatment Condition",
       x = "Prior Knowledge",
       y = "Frequency") +
  facet_wrap(~ treat) +
  theme_bw()




